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INTRODUCTION 

This  report  describes  the  continuation  of  two  efforts.  The  first  is  a  computing  scheme  de¬ 
veloped  by  the  authors  for  the  analysis  of  factorial  designs.  An  IBM  704  Electronic  Data  Proc¬ 
essing  Machine  program  using  this  scheme  was  submitted  to  the  SHARE  Society  in  1958,  and  a 
short  paper  discussing  this  scheme  was  presented  to  the  14th  national  meeting  of  the  Association 
for  Computing  Machinery  (Ref.  1). 

The  second  effort  is  an  attempt  to  promote  a  general  linear  approach  to  numerical  and  sta¬ 
tistical  analysis.  A  report  describing  this  approach  in  the  smoothing  of  experimental  data  has 
already  been  published  (Ref.  2). 

Three  concepts  form  the  basis  of  these  efforts:  the  Gram-Schmidt  construction,  the  linear 
method  of  least  squares,  and  the  application  of  orthogonal  polynomials.  References  3.  4,  and  5, 
and  many  other  texts,  also  introduce  these  well-known  concepts. 

The  notation  of  vectors,  matrices,  and  tensors  is  used  in  place  of  the  summation  notation 
frequently  (but  aot  universally)  used  in  textbooks.  It  is  believed  that  the  notation  used  is  more 
practical  at  the  working  level,  with  the  advent  of  FORTRAN  and  ALGOL,  and  it  certainly  allows 
a  much  shorter  text. 

Computer  programs  for  factorial  and  statistical  analysis  are  presented  in  Appendixes  A.  B, 
and  C:  they  demonstrate  a  method  rather  than  finished  programs.  They  are  finished  to  the  extent 
that  results  are  summarized,  but  there  are  many  wavs  of  summarizing,  depending  principally  on 
the  desires  of  the  user.  Because  of  the  versatility  of  the  FORTRAN  language,  the  prospective 
user  is  invited  to  summarize  results  in  the  manner  he  sees  fit. 

This  work  is  presented  principally  to  those  who  are  comput-T-oriented.  For  this  reason, 
there  is  a  great  deal  of  elaboration  on  some  of  the  basic  concepts  and  philosophy  of  statistics. 

On  the  other  hand,  it  is  assumed  that  the  reader  is  familiar  with  the  algebra  of  vectors  and 
mat-ices  and  their  geometric  interpretations. 


VECTORS,  MATRICES,  AND  THE 
GRAM-SCHMIDT  CONSTRUCTION 

Linear  hypothesis  statistics,  which  includes  regression  and  factorial  analyses,  may  be  de¬ 
veloped  by  a  method  (the  ortliogonalization  of  matri'-es).  a  theorem  (the  Pythagorean),  and  a 
philosophy  (inference).  This  development  is  best  made  with  a  uniform  notation.  Therefore,  the 
following  definitions  are  made. 

Matrices  will  be  denoted  by  capital  letters  such  as  V.  )  ,  .1,  and  D.  und  scalar  quantities 
will  be  presented  by  small  letters  such  as  a,  b.  and  c. 

The  elements  of  a  matrix  will  be  represented  by  subscripted  small  letters  such  as  x-j  and 
v .. 

•  i 

The  ith  row  of  a  matrix  A  will  be  denoted  bv  x i  ,  and  the  /th  column  of  a  matrix  V  will  be 
denoted  by  t  ..  A  row  or  a  column  (or  anv  singly  subscripted  variable)  may  be  considered  to  be 
a  vector  in  space,  the  elements  of  the  row  or  column  being  the  ordered  components  of  the  vector. 
This  is  the  geometric  viewpoint. 
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Two  columns  of  a  matrix  are  said  to  be  orthogonal  if 
x.j  *  x.k  "  0  /  ^ k 

where  the  vector  dot  product  notation  above  is  defined  thus 
n 

x.j  ‘  *.*  *£(*.7)(xi*) 


(i) 


(2) 


i-  1 


From  the  geometric  viewpoint,  the  vectors  x  •  and  x  >  meet  at  right  angles  in  n-dimensional 
space.  ■' 

An  individual  column  of  a  matrix  ,Y  is  said  to  be  normalized  if 


An  unnormalized  column  (provided  it  contains  at  .east  one  finite  entry)  may  be  easily 
normalized  by  dividing  each  element  of  the  column  by  the  length  of  the  column  vector  as  shown: 

1  <  i  <  n 


x'.. 
1 1 


yi* 


and  where  the  length  is  defined  as  follows: 


i-  1 


(4) 


From  the  geometric  viewpoint,  u  normalized  column  is  a  vector  of  unit  length.  The  vector  need 
not  be  oriented  along  any  axis  but  may  have  components  along  all  axes.  Fquation  f  rewritten 


-  x 


X. 

•I 


i  I 


is  a  statement  of  the  Pythagorean  theorem  generalized  to  m  dimensions.  The  square  of  the 
hypotenuse  (vector  length)  is  equal  to  the  sum  of  the  squares  of  the  other  sides  (components). 


A  matrix  ,\  is  said  to  be  orthonormal  (with  respect  to  columns)  if  each  column  is  normalized 
and  if  each  column  is  orthogonal  to  every  other  column.  From  a  basic  theorem  of  linear  algebra, 
this  is  possible  only  if  the  number  of  columns  does  not  exceed  the  number  of  rows  (Ref.  3). 


An  important  property  of  any  orthonormal  matrix  is  that  premultiplication  by  its  own  transpose* 
yields  an  identity  matrix. 


P'  P  I 


(5) 


This  may  Se  demonstrated  ns  follows: 


- 1 

iF 

1  0 

l-fl  0  fl 

f  0 

m  m 

o  1 

Note  that  it  is  the  second  matrix  which  is  orthonormal  with  respect  to  columns. 


I 


I  h<‘  same  matrix  isilh  rows  and  column* 


inlmhanpod 
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Conversely,  the  product  P  P'does  not  generally  yield  an  identity  matiix  unless  P  is  square. 
This  also  may  be  demonstrated 


- _ 1 

- 1 

O' pi 

1 

3 

i ' 

"5 

J1  0 

a 

o  .Jl 

1 

1 

1 

’  3 

\2 

N  2 

3 

3 

3 

1 

6 

m 

I 

3 

s 

i .}  n  2 

3 

0 

A  matrix  .V  may  be  transformed  into  an  orthonomial  matrix  P  by  the  Gram-Schmidt  construc¬ 
tion  (Ref.  3).*  The  following  paragraphs  describe  this  construction.  The  matrix  is  first  trans¬ 
formed  into  a  matrix  having  the  orthogonal  propertv 


v  — ►  V 

Thence  it  is  transformed  into  a  matrix  which  also  has  the  normal  property 

V  —  P 

The  first  column  of  the  orthogonal  0  matrix  is  formed  by  simply  using  the  first  column  of 
the  V  matrix'1 

q  ,  -  x  ,  (6) 

The  second  column  of  the  Q  matrix  is  formed  by  taking  the  second  column  of  the  A  matrix 
and  adding  to  it  that  portion  .  v  ,2  .  of  the  first  column  of  the  Q  matrix  which  makes  the  new 
column  orthogonal  to  the  first  column  of  the  Q  matrix.  The  following  vector  equation  expresses 
this  required  condition 


<7.  I  ’  (*.2  +  r  12  ?.1)  ;  0 

The  solution  yields 

f  "*.»  •  Ll 

12  7 .  i.  •  <7,  i 

With  this  result,  the  second  column  of  the  Q  matrix  is 

f.2  ~  *.2  +  r  1 2  1  ^ 

Tbe  third  column  of  the  (J  matrix  is  formed  by  taking  x  3  and  adding  portions  cl3  and  c23 
of  the  q  ,  and  q  2  columns  which  make  q  3 orthogonal  to  both  q  ,  and  q  2-  Two  simultaneous  vector 
equations  express  these  desired  conditions 


<7.1 

-  (\3  +  rt3<?. 

1  4  c  2  3  *7.2^ 

0 

(9) 

<7.2 

*  (t.3  +  r13  <7. 

1  +  r23  *7.2^ 

0 

(10) 

2 

‘‘Provided  the  columns  have  the  property  of  linear  independence  A  further  elaboration  of  this  is  given 
after  the  mechanics  of  the  construction  is  presented 

1  Remembering  that  the  dot  notation  is  defined  thus 
*■  l  ~  *  1 1  I  <  *  <  " 
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Since  q  t  •  q  2  -  0..  c  ,  3  and  c  23  are  quickly  evaluated  to  be 
-9,1  *  *.3 

mi 

9.2  *  9.2 

With  these  results,  the  third  column  of  the  orthogonal  matrix  is 

?  }  =  *  3*  f13?  J  1  C  23  9  2  03) 

The  Gram-Schmidt  construction  may  be  easily  generalized  so  that  the  expression  for  the  /th 
column,  q  ■.  may  be  written 
/  I 

,r  !</<«  (M) 


I  <  i  < 


This  may  be  expressed  by  the  matrix  relation 

^  -  ,Y  *  QC 


where 


.  =  o 
*/ 

/<* 

~9.i  *  \i 

ki  9.*  *  H.h 

/>* 

I  £  /  S  m 
I  £  A  <  m 


and  where  m  is  the  number  of  columns  and  n  is  the  number  of  rows  of  the  A  and  Q  matrices. 

The  orthogonal  matrix  Q  mav  now  be  transformed  into  the  orthonormal  matrix  P  by  dividing 
each  element  q-j  by  the  length  of  the  column  vector  to  which  it  belongs 

Pii  ^9./ 1  '  'lij  *  <i<* 

I  <  j  <  m 

This  too  may  be  expressed  by  a  matrix  relation 


where 


P 

QR 

0* 

--  0 

i*k  ) 

1  <  /  <  m 

r/* 

/  *  I 

1  <i  <  m 

Normalization  may  not  be  possible.  This  is  obviously  the  case  where  all  elements  of  a 
column  are  zero.  When  such  a  column  vector  appears  in  a  matrix  during  Gram-Schmidt  construc¬ 
tion,  the  original  matrix  lacks  the  required  property  of  linear  independence.  In  fact,  this  is  an 
excellent  practical  test  of  linear  independence. 

Two  important  relations  have  been  developed  above 

Q  Y  *  QC  OS) 

P  QR  06) 

The  first  mav  be  rearranged  as  follows: 

Qil  ~  C)  -  ,Y  07) 
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(/  -  C)  is  a  triangular  matrix  with  diagonal  terms  of  unity;  its  determinant  is  nonzero  (in  fact, 
unity)  and.  therefore,  its  inverse  exists,  and  Eq.  17  may  be  rewritten 

Q  .  X  (/  -  O'  *  (18) 

From  Fq.  16.  P  is  found  to  be 

P  -  V  (/  -  cr 1  H  (19) 

or  P  -  ,V  D 


where  l)  -  (/  -  C)~  1  R 


Equation  19  demonstrates  that  the  Cram-Schniidt  construction  is  equivalent  to  postmultiplication 
of  the  original  matrix  bv  another  matrix.  I). 

There  is  also  a  geometric  viewpoint  of  the  construction.  The  columns  of  the  A  matrix  repre¬ 
sent  vectors  of  \ar\ing  lengths,  meeting  each  other  at  various  angles  in  n-diinensional  space. 

Each  vector  is  considered  in  the  arbitrary  order  in  which  they  appear  in  the  A  matrix,  and  its 
components  which  lie  along  preciously  chosen  vectors  are  removed.  The  first  vector  obviously  has 
no  other  vector  to  be  compared  with  and  is  left  alone. 


This  can  be  illustrated  graphically  in  the  two-dimensional  case  (Fig.  1).  Consider  two 
vectors,  x  j  and  t  2>  described  by  a  matrix  X.  where  x  is  the  component  of  x  j  on  the  it  axis, 

x,  axis 


H(,  |  Geometric  Representation  of  the  \  anil  Q  Matrices  in  Two  Dimensions,  Showing 

the  (>rani-S<  hmiift  Construction. 
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and  x 21  ia  the  component  of  x  ,  on  the  x2  axis,  x,2  is  the  component  of  x  2  on  *l*e  *1  *■***.  and 
*22  is  the  component  of  x  2  on  the  x2  axis.  In  the  Gram-Schmidt  construction  the  column  x.j  is 
accepted,  and  its  elements  are  put  into  the  first  column  of  the  matrix  Q. 

Then  the  vector  x  2  is  taken,  and  its  component  (cxt()  along  .v,j  is  subtracted  to  form  the  vector 
q  2>  The  components  o(q2  along  the  x(  and  *2  axes  now  form  the  second  column  of  an  orthog¬ 
onal  matrix  Q.  The  normalizing  process  (not  shown  in  Kig.  1)  shrinks  q  j  and  q  2  to  unit  lengths 
but  preserves  their  respective  directions,  which  now  form  a  right  angle  with  each  other. 

On  digital  computers,  the  matrices  P  and  D  are  determined  from  X  as  follows:  The  P  matrix 
is  formed  in  the  same  memory  area  as  the  matrix.  Kach  succeeding  column  vector  of  the  X 
matrix  (except  the  first)  is  orthogonalizcd  separately  with  respect  to  each  succeeding  column  in 
the  manner  previously  described.  The  s.  parate  columns  of  the  matrix  may  then  be  normalized. 
Meanwhile,  the  same  operations  performed  on  X  are  also  performed  in  another  memory  area  orig¬ 
inally  containing  an  identity  matrix  /.  It  is  obvious  that  as  .V  is  transformed  into  P ,  I  is  trans¬ 
formed  into  D,  thus 

X  — ►  P 

i  — >  n 

because,  as  noted  above,  this  construction  is  equivalent  to  postmultiplication  by  the  matrix  D. 

A  numerical  example  of  these  transformations  appears  in  Appendix  D. 


REGRESSION  ANALYSIS  AND  THE  LINEAR 
METHOD  OF  LEAST  SQUARES 


The  linear  method  of  least  squares  is  closely  related  to  the  solution  of  simultaneous  linear 
equations  by  the  method  of  orthogonal i/ation  which  will  be  briefly  described. 

The  following  system  of  linear  equations, 


al  x\  I  +  “2*12  1 


a 1  *21  +  a2r22 


•  *+  anXln  "7, 


4  ax. 
n  2/i 


°1 +  fl2*„2+ 


tax  -  y 
n  nn  J  n 


may  be  represented  by  the  matrix  equation 


where  the  .V  and  1’  matrices  are  known,  and  .1  is  to  be  determined. 

With  the  results  of  the  first  section,  a  solution  may  be  found  as  follows: 
PD-'A-.Y  (from  Kq.  19) 

P'  P  /T1  I  -  P'  ) 

[)-'  A  -P"  > 


6 


(from  Kq.  5) 


(22) 

(23) 

(24) 
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hence 


A  =  DP'Y  (25) 

A  FORTRAN  coded  subroutine  for  simultaneous  linear  equations  using  this  method  appears 
in  Ref.  2. 

At  this  point,  it  is  useful  to  look  at  the  geometric  viewpoint  of  the  system  of  linear  equations 
and  its  solution  as  outlined  above.  1 ,  of  course,  should  be  considered  a  vector  in  n-dimensional 
Euclidean  space,  with  a  direction,  a  length,  and  components  y-.  As  before,  the  dot  product  of 
)’  on  )’  should  be  recognized  as  an  expression  of  the  Pythagorean  theorem: 

>■  .  Y  ,  |  >i  a  =£y.2  (26) 

where  j  ) '  |  is  the  length  of  the  vector.  The  y^.  are  components  of  )  when  referred  to  a  particular 
referenc  e  frame,  namely,  that  defined  by  the  identity  matrix  /  of  order  n. 

Note  that  Eq.  25  mav  be  rewritten  R  -  Pr  ) .  where  R  •-  Since  P  is  orthonormal,  the 

elements  of  R  are  also  components  cf  )7*  This  time  the  components  h(.  are  along  unit  vectors 
whose  components,  in  tum.  are  found  in  the  columns  of  P,  and  hence  there  is  a  new  reference 
frame  P  .  The  Pythagorean  theorem  may  also  be  stated  in  terms  of  this  reference  frame 

l>  I2  £V2  -  l«l2  (27) 

Inherent  in  this  development  is  the  fact  that  an  orthonormal  matrix  is  a  vector  basis ,  or 
Cartesian  reference  frame,  in  which  the  )  vector  retains  both  its  magnitude  and  direction  In 
other  words,  the  vector  R  in  the  P  reference  frame  is  identical  with  the  )  vector  in  the  /  refer¬ 
ence  frame. 

The  solution  process  can  be  described  in  this  light.  Consider  the  matrix  equation  A  A  ■*  Y. 

•I  may  be  called  a  representation  of  I  .  because  it  may  be  determined  from  )'.  Even  more,  it 
represents  }'  in  a  manner  very  pertinent  to  .V.  so  it  may  be  stated  further  that  A  is  a  representa¬ 
tion  of  )  in  the  A  space.  Hut  it  is  u  distorted  representation  of  )',  because  the  A  space  is 
generallv  distorted.6  If  this  distortion  is  removed  by  Gram-Schmidt  construction,  there  is  ob¬ 
tained  a  space,  P.  in  which  )  has  an  undistorted,  easily  calculated,  understandable,  representa¬ 
tion,  R.  But  R  is  not  what  is  wanted,  so.  to  get  I.  one  must  reintroduce  the  distortion.  This  is 
done  in  the  final  step  of  the  solution  I  -  DR  (from  Eq.  25).  Thus  the  term  “distortion,"  used 
very  loosely  above,  is,  nevertheless,  precisely  and  quantitatively  expressed  by  the  matrix  D. 

To  turn  now  from  the  geometric  viewpoint,  consider  the  case  where  the  columns  of  the  X 
matrix  do  not  possess  linear  independence.  This  is  revealed  when  a  column  of  zeros  appears 
during  the  Gram-Schmidt  construction.  When  this  happens,  the  vector  basis  is  incomplete,  be¬ 
cause  there  are  fewer  vectors  than  dimensions.  The  vector  basis  is  also  incomplete,  for  the 
same  reason,  when  there  are  fewer  columns  than  rows  in  A.  This  will  be  considered  to  be  the 
general  case  of  the  incomplete  vector  basis,  because  zero  columns  which  appear  during  Gram- 
Schmidt  construction  may  be  simplv  deleted  with  no  loss  of  generality. 

It  will  first  be  shown  why  this  situation  (fewer  columns  than  rows)  actually  occurs.  Suppose 
one  is  conducting  a  scries  of  experiments  to  measure  the  value  of  y.  the  dependent  variable. 

I'nder  experimental  control  are  m  independent  variables  x  ■ .  By  physical  law  or  statistical  hy¬ 
pothesis,  y  may  vary  linearly  with  the  x .  as  follows: 

r=al*l  .  a2*2  ♦  .  .  .«■  amx„  (28) 

Equation  28  is  called  the  “regression  model  ." 

^Reference  .1,  or  another  basic  text,  may  he  consulted  to  prove  this  point 

It  is  simply  two  different  methods  of  locating  its  coordinate  points 

6  Not  orthonormul 
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Obviously,  one  must  perform  at  least  m  experiments  to  determine  the  a,.  However,  one  may 
proceed  further  and  perform  n  experiments,  where  n  is  any  number  greater  than  m.  Attaching  the 
subscript  i  for  the  ith  experiment,  one  gets  the  following  set  of  “overdetermined"  equations 

a.  *il  *  a2Xi2  +  •  *  •  +  am  Xim  "  7i  ^  «  <*> 

which  may  also  be  represented  by  the  matrix  equation 

XA  ,  Y  (21) 

noting  that  ,¥  now  contains  more  rows  than  columns. 

While  it  is  not  surprising  that  the  experimenter  did  not  have  to  stop  after  m  experiments,  it 
may  be  surprising  that  the  overdetermined' equations  above  may  be  solved  in  the  same  manner 
as  outlined  in  Eq.  21  through  25.  A  careful  examination  will  show  that  no  step  is  invalid. 

The  preceding  is  contrary  to  usual  experience  in  solving  systems  of  linear  equations.  If 
the  X  matrix  is  not  linearly  independent,  no  solution  is  possible  by  the  usual  methods.  But 
here  there  is  a  solution  where  the  A’  matrix  is  not  even  square,  let  alone  linearly'independent  . 
Before  examining  the  meaning  of  such  a  solution,  a  numerical  illustration  will  be  given  . 


I,et  the  dependent  variable  be  the  yield  of  the  principal  product  in  a  chemical  experiment. 
One  is  interested  in  how  the  yield  varies  with  pressure  (x2),  temperatrre  (x$),  and  amount  of 
catalyst  (xj).  Six  experiments  provide  the  following  data: 


Experiment  No. 

Pre 

ssure 

Temperature 

Catulyst 

Yield 

U) 

( 

*i2> 

U.V 

1 

1 

300 

.001 

108 

2 

1 

too 

.001 

111 

3 

1 

600 

002 

120 

1 

2 

350 

.001 

99 

5 

2 

5.50 

.002 

115 

6 

2 

650 

.002 

117 

The 

following  regression 

model  can  then  be 

set  up: 

al 

*«l  ’  U2Xi2  ' 

‘VV't  f 

a .  x. 
4 

\-y 

i  1  2 

i  < 

5 

(30) 

In  this 

rase 

(which  is  very  common),  a  dummy  X|  vari 

able  must 

be  defined 

which  takes 

the 

value 

unity  for  all 

values  of 

The  data 

and  model  are  now 

expresser 

in  the  foil 

ow  ing  matrix 

: 

*  i 

1  300 

001’ 

■»  ^ 

a, 

'108' 

i 

1  100 

.001 

a2 

Ill 

i 

1  600 

1)02 

U3 

120 

(31) 

i 

2  350 

.001 

“j 

99 

i 

2  .550 

.002 

L  - 

115 

.  i 

2  650 

1)02. 

.117. 

P  ar 

d  I)  mat 

rices  are  now 

found  by 

Gram-Schmidt  construction,  and  the  .1  matrix  is  found  bv 

Kq. 

25  to  be 

(99.H7.  -7.2 

.50.  .025. 

7.500) . 

What  is 

the  1  vector 

in  this  case.’' 

The 

.1  vector 

s  an  incomplete  representation  of  the 

}’ 

vector  in  the 

A  space.  It 

has  only 

four 

■omponents:  to 

be  compl 

etc,  it  shou 

Id  have  six 

8 

Further 

insight  may 

ie  gainec 

by  re 

turn 

ing  to  the 

geometric 

viewpoint 

and  to  the  soluti 

in 

proc 

ess 

itse 

If.  The  Gram 

-Schmidt 

const 

ruct 

on  has  de 

veloped  a 

matrix  P.  which  contai 

is  four 

"  M 

>rc  row*,  than  column*  in  the 

\  matrix;  «* 

>|uix.il<ntly 

mon* 

lat  a 

than  paramrt 

rrs  in  an  exp 

eriment. 

Hln 

Mill 

*tii  al  termino 

ojey.  «  omponc*nt> 

are 

(re. 

ilom . 
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columns  in  the  example.  These  represent  four  vectors  meeting  at  right  angles  in  six-dimensional 
space.  Therefore,  when  the  solution  step  R  -  P ')  is  performed,  only  four  of  the  six  possible 
components,  b ■,  are  computed.  It  would  be  possible  to  compute  two  more  components  by  arbi¬ 
trarily  adding  two  more  orthogonal  vectors  to  the  P  space.  However,  these  two  new  vectors 
would  have  no  definite  relation  to  the  independent  variables  of  the  experiment,  so  this  addition 
would  not  be  justified.  One  must  accept  the  actual  situation  . 

R,  less  two  components,  while  not  a  complete  representation  of  )'  in  the  P  space,  is  a  com¬ 
plete  representation  of  some  other  vector,  f .  which  does  not  have  components  along  the  two 
missing  P  vectors,  f  can  be  found  by  the  reverse  transformation  f  =  PR,  which  transforms  R 
back  onto  the  /  reference  frame.  (  This  expression  mav  be  shown  to  be  equivalent  to  f  =  XA, 
where  A  is  obtained  from  Kq.  2.3.) 

In  the  example.  V  and  Y  compare  os  follows: 

f  > 


107.2 

103 

109.7 

ill 

122.2 

T20 

101.2 

99 

113.7 

11.3 

116.2 

117 

It  is  not  an  exact  solution,  but  it  is  close,  and  there  is  good  reason  why  it  should  be  close.  If 
the  physical  law  is  correct  and  there  are  no  experimental  errors,  one  would  expect  the  Y  vector 
to  have  onlv  four  components  in  its  A  representation,  namely,  those  of  the  regression  model. 

The  manifestation  of  other  components,  as  seen  in  the  difference  between  )  and  }",  would  be  ex¬ 
pected  only  if  there  were  experimental  error  or  if  the  model  were  incorrect.  Thus,  it  is  not  only 
possible  to  have  an  overdetermined  set  of  equations,  but  it  can  be  very  desirable,  because  it 
can  elucidate  sources  of  error.  It  will  now  be  shown  that  this  is  not  only  a  solution  to  the  over¬ 
determined  set  of  linear  equations,  but  that  it  is  the  best  solution  in  the  sense  of  least  squares. 

The  problem  of  least  squares  can  be  stated  from  the  geometric  viewpoint: 

The  best  rep--esentat  ionQ  of  a  vector  )  of  dimension  n,  in  a  subspace  P  of  dimension  m, 
smaller  than  n,  is  found  bv  constructing  the  vector  J  ,  having  components  only  in  P,  such  that 
the  square  of  the  distance  between  }  and  )  is  a  minimum. 

The  distance  between  two  vectors  is  the  length  of  the  vector  difference,  denote  J  }’  -  f  J  . 

The  evaluation  of  this  quantity  follows.  The  unknown  vector  )  has  components  Jj,  i . 

in  the  P  space,  which  are  chosen  to  fulfill  the  criterion  above.  If  additional  vectors  are  added 
to  /Min  practice,  a  trial  .md  error  process)  so  as  to  form  a  complete  orthonormal  vector  basis, 

the  components  of  the  known  vector  }  may  be  exactly  determined,'11  bt,  b i . bm,  t . 

bn,  where  the  last  n  -  m  components  do  not  belong  to  the  original  incomplete  space.  The  compo¬ 
nents  of  the  vector  difference  I  —  f  are  found  bv  the  rule  of  vector  subtraction 


<bx~b,).(b2-bS . (bm  -  bj.bmtl . bn 

The  square  of  the  vector  difference  (which  is  to  be  minimized)  is  found  by  the  theorem  of 
Pythagoras 

-  c*.  -  6,>=  •...<»„■■  t.)’  •  «2) 


9  In  the  sense  of  least  squares. 
As  before,  R  =  P'  I . 
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The  minimum  occurs  when  each  h,  is  set  equal  to  b ,•  for  all  i  up  to  m;  this  can  be  seen  by  ob¬ 
serving  that  each  of  the  leading  terms  has  a  minimum  possible  value  of  zero.  Thus,  the  criterion 
of  least  squares  is  fulfilled  when  )'  consists  of  the  components  of  Y  in  the  incomplete  P  space. 
This  is  exactly  what  the  solution  process  above  does,  and  it  will  henceforth  be  referred  to  as 
the  linear  method  of  least  squares. 

The  foregoing  is  not  the  classical  approach  to  least  squares,  nor  does  it  yield  the  more  fa¬ 
miliar  least  squares  equations.  It  is  a  curiosity  in  statistics  that,  while  much  of  the  theory  is 
based  on  linear  hypotheses,  much  of  the  analysis  is  presented  in  a  cumbersome  quadratic  form  . 
The  subject  of  least  squares  is  easily  presented  with  the  aid  of  the  Pythagorean  theorem,  ns 
shown  above.  In  the  classical  approach,  differential  calculus  and  summation  algebra  are  used. 
To  arrive  at  the  quadratic  form  of  the  least  squures  equations,  modify  the  linear  equation  Ad  * 
by  premultiplying  both  sides  by  .V/,  which  yields  (A  ’  AVI  -  (A'/)'),  which,  expanded,  yields 


This  is  the  form  presented  in  many  texts. 

A  sample  coding  for  regression  analysis,  using  the  linear  method,  appears  in  Appendix  A. 
Statisticians  use  the  inverse  matrix  elements  from  the  quadratic  form  to  perform  t  tests,  which 
show  whether  the  data  conform  to  the  regression  model  or  not.  The  reader  is  referred  to  stand¬ 
ard  works  (  such  as  lief,  t)  for  a  description  of  these  tests.  I’he  statistician  is  not  deprived 
of  these  elements  in  the  linear  method,  for  one  can  show  that  the  inverse  matrix  (A*  V) "  is 
equal  to  the  product  DP'  from  the  linear  method.  The  sample  coding  in  Appendix  A  includes 
this  calculation. 


ORTHOGONAL  POLYNOMIALS  AND 
CURVILINEAR  REGRESSION 

In  the  previous  secti  r-s  it  was  seotvn  that  the  square  matrix  A  of  a  matrix  equation  may  be 
transformed  into  a  matrix  t'  which  has  highly  desirable  algebraic  and  more  reudilv  understood 
geometric  properties.  The  column  vectors  of  the  P  matrix  are  of  unit  length  and  are  mutually 
perpendicu lar  to  one  a  i  -'i.er.  •  urthe-more,  the  set  of  vectors  may  be  incomplete,  in  which  case 
there  are  fewer  vectors  item  data  points. 

The  elements  of  the  or.hogt  nai  vectors  are  called  by  statisticians  "orthogonal  polynomials, 
for  the  vector  elements  are  actually  values  taken  by  the  set  of  polynomial  functions 
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Pi\  *  'll*,! 

Pn-'.,*,,*'  »*«  !  'i-'S”  ls" 

P<  3  *  ^1  3  *|  I  +  ^23*»2  +  ^33  *13 
These  equations  are  an  expansion  of  the  matrix  product,  P~XD,  from  the  preceding  section. 

The  most  common  application  of  orthogonal  polynomials  occurs  in  curvilinear  regression . 
Curvilinear  regression  is  a  special  case  of  regression  analysis  based  on  the  following  model 


•  •  o. 


xmf{  I  <  i<  n 


=  +a2*>t03  V - m-i 

Again,  this  may  be  represented  by  a  matrix  equation 


VI 

> 

where  the 

i/th  el 

and  where  the  singly 
/'-I  is  an  exponent.  : 
A  as  follows: 

Pi  1 

Pi  2 

d\2  * 

Pi  3  - 

d  1  3  * 

(35) 


(21) 


(36) 


d22  Xi 


(37) 


In  the  functions  above.  p(  )  is  a  constunt.  p-2  is  a  linear  function  of  r(,  p(  ^  is  a  quadratic 
function  of  and  so  on.  A  very  important  facet  of  these  equations  should  be  noted  here. 
Although  the  members  on  the  right  are  poly  nominal  functions  (of  linear,  quadratic,  and  higher 
degree  varieties),  the  dependent  left-hand  members  (being  only  the  set  of  func  tional  values)  may 
bv  viewed  as  vectors  and  represented  in  a  linear  vector  space.  This  will  be  an  aid  to  under¬ 
standing  how  one  can  speak  of  higher  order  variations  with  a  lineur  theory  . 

Continuing  to  the  next  step  of  the  linear  method 

li  P'Y  (38) 

This  equation  may  be  broken  down  as  follows: 


r.i-r. 
V  p.2  *r. 


(39) 


^  3  p.  .3  "  y. 


The  significance  of  the  />..  as  the  components  of  )  in  the  P  reference  frame,  has  already  been 
discussed.  In  addition  to  the  algebraic  and  geometric  meaning,  the  b  ■  will  be  shown  to  take  on  an¬ 
other  meaning.  It  is  about  the  h.  components  of  }  that  the  statistician  begins  making  inferences. 
It  is  inference  that  separates  the  statistics  from  the  physics  and  mathematics  of  the  experiment 
and  of  its  model.  Kquation  39  will  serve  as  an  illustration.  I’be  term  b2  may  be  inferred  to  be  the 
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variation  of  yi  due  to  a  linear  variation  of  x.;  b3  may  be  inferred  to  be  due  to  the  quadratic  vari¬ 
ation,  and  ao  on  through  higher  ordera. 

That  inference  ia  not  atrictly  mathematical  may  be  aeen  by  the  word  “due.”  There  ia  almost 
a  cauaal  implication  here,  although  statistics  are  not  supposed  to  show  causesM  What  is  meant 
is  that  the  bj  (which,  mathematically,  are  only  linear  vector  components)  are  inferred  to  show  a 
true  functional  relation  (not  necessarily  linear)  between  the  dependent  and  independent  variables. 
Under  certain  frequent  conditions,12  it  is  possible  to  show  that  the  b .  are  indeed  the  best  esti¬ 
mates  of  the  true  functional  relations,  in  which  case  the  statistician  has  great  confidence  in  his 
inferences. 

Continuing  with  the  linear  method  for  the  case  at  hand  (curvilinear),  the  last  step  is  to  trans¬ 
form  the  B  vector  into  the  A  vector,  by  premultiplying  with  the  D  matrix,  and  thus  complete  the 
determination  of  the  regression  coefficients.  This  last  step  may  be  considered  to  be  an  anti¬ 
climax.  Before  the  last  transformation,  there  existed  an  ideal  Euclidean*"*  space,  where  the  bj 
were  merely  components  of  Y  in  the  P  reference  frame.  The  last  transformation  not  only  destroys 
this  situation  but  also  destroys  the  inferences  as  well,  because  the  magnitude  of  the  a(.  depend 
not  only  on  the  linear,  quadratic,  cubic,  etc.,  variations  of  y  with  x,  but  also  on  the  amount  of 
distortion  in  the  X  matrix. 

If  the  last  transformation  is  deferred,  there  are  not  only  components  about  which  inferences 
can  be  made,  but  there  is  also  a  simple  device  at  hand  to  test  the  validity  of  the  model.  In  the 
preceding  section,  it  was  shown  that  Y  may  have  components  which  do  not  appear  in  the  P 
reference  frame.  The  hypothesis  may  be  made  that  these  components  are  experimental  error.  If 
the  error  is  small,  the  components  of  Y  about  which  functional  relations  of  y  with  x  were  in¬ 
ferred  should  be  larger  than  those  about  which  functional  relations  were  not  inferred.  The 
following  comparison  function  may  therefore  be  set  up 


k~m+  1 


(40) 


which  is  called  an  F  ratio.  Note  the  b ^  in  the  denominator  are  those  components  of  y  which  do 
not  appear  in  the  P  space. 

An  F  distribution  is  a  function  of  the  number  of  inference  components  in  the  numerator  and  the 
number  of  error  components  in  the  denominator.  When  an  F  is  larger  than  a  preset  percentage 
(say  95°S)  of  all  other  possible  F s  that  are  determined  by  the  distribution,  the  statistician  says 
that  a  true  functional  relation  exists.  In  normal  practice,  the  computer  calculates  the  F  ratios, 
and  the  statistician  draws  inferences  and  performs  the  F  tests.  F  and  t  tests  are  part  of  the 
probability  aspects  of  statistics,  and  it  is  not  in  the  scope  of  this  report  to  pursue  them  further. 
The  reader  is  referred  to  Ref.  4  and  other  basic  works  on  statistics. 


**  For  an  interesting  discussion  of  this,  see  How  To  Lie  With  Statistics,  by  D.  Hoff,  W.  W.  Norton  and  Co., 
New  York,  1954.  142  pp. 

*2  Namely,  notmal  (Gaussian)  distribution  of  error. 

‘■’That  is,  undistorted,  orthonormal. 
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The  computational  aspects  of  curvilinear  regression  will  now  be  considered.  It  will  be  shown 
that  the  calculations  of  the  linear  method  can  be  performed  stepwise.  As  the  name  implies,  step¬ 
wise  calculation  allows  the  regression  coefficients  to  be  calculated  and  tested  one  at  a  time. 

Such  a  process  is  desirable  because  the  curvilinear  regression  model  approximates  a  Taylor  s 
expansion  of  y  on  x 

yr ■■  c  ,  +  o2  x  f  «,  x2  ♦  .  .  .  ♦  a,  x"~  1  <41) 

where 


LW  J«  «•••»■ 


(42) 


Because  the  reciprocal  factorial  term  in  «■  rapidly  approaches  exceedingly  small  numbers,  it  is 
expected  that  at  some  term  in  the  model,  the  remaining  functional  relation  between  t  and  y  is 
small  compared  to  the  experimental  error.  4  stepwise  solution  can  locate  this  term  and  can, 
furthermore,  yield  regression  coefficients  for  all  intermediate  degree  models  without  numerical 
repetition. 


Having  established  the  motive  for  using  stepwise  regression,  it  will  be  shown  how  it  works 
for  digital  computers.  One  possible  method  is  to  use  published  values  of  orthogonal  polynomials, 
but  this  is  not  usually  practical  on  digital  computers.  Only  the  orthogonal  polynomials  for  equally 
spaced  intervals  of  the  independent  variable  are  usually  published.  ,Also,  the  D  matrix  is  very 
hard  to  determine,  requiring  almost  as  much  computation  as  the  Grain-Sehmidt  construction.  A 
more  practical  approach  to  the  use  of  digital  computers  will  be  described. 

Stepwise  regression  is  the  technique  of  performing  all  steps  of  the  linear  method  by  analyzing 
one  column  of  the  A  matrix  at  a  time,  according  to  the  following  outline,  which  describes  the  steps 
for  the  y t h  column 


1.  Korm  the  yt h  column  of  the  A  matrix 


,t . .  -  \/~  1  1  <  i  <  n  (43) 

2.  Orthogonal ize  this  column  ^with  respect  to  previous  columns  by  the  Gram-Schmidt  construc¬ 
tion.  and  normalize  to  find  the  yth  column  of  the  P  and  P  matrices 


Pip  tlij  1 

3.  bind  the  yth  component  of  )  in  the  P  space'4 

(44) 

h .  p  ’  v 
/  W  • 

(45) 

t.  Compute  the  current  F  ratio 

h)2 

F  11.  n-j)  - l - 

(46) 

A  j  *  1 

3.  The  F  ratio  of  step  t  may  be  used  to  test  for  significance  of  the  yth  functional  relation. 
At  any  time  before  column  m  is  reached,  the  statistician  may  decide  to  stop  the  stepwise 
regression  on  the  basis  of  this  test,  or  perhaps  on  the  basis  of  the  residual  standard 
deviation  discussed  below.  It  is  difficult  to  get  agreement  on  a  proper  criterion. 


It 


I'he  vector  v . 


is  identical  I  o  the  )  matrix 
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6.  At  any  time,  the  conventional  regression  coefficients  for  a  /'-I  degree  fit  can  be  found 
from  the  current  states  of  the  D  and  li  matrices 

°i  dij  •  bj  (47) 

This  technique  is  used  in  the  sample  program  shown  in  Appendix  B. 

The  denominator  in  Kq.  16  is  a  useful  statistic  called  the  “variance,*  s2.  of  the  regression 


n 

1  V 

«-/  /  ^ 


(48) 


k j¥\ 


The  square  root  of  this  quantity,  namely,  the  residual  standard  deviation,  is  an  estimate  of  the 
experimental  error. 

Note  that  Kq.  46  and  48  contain  components  of  li  which  have  yet  to  be  computed.  It  is  known, 
n  n 

however.  that^T^hy2  ^™'vj-2.  from  the  Pythagorean  theorem,  and  therefore  the  following 
r*  ,  i 

relation  may  be  used  to  compute  the  required  components: 
n  n  I 


2>!  2}<! 

TjT i  k  i  *7 


49) 


THE  DESIGN  OF  EXPERIMENTS 

In  the  preceding  section,  it  was  noted  that  useful  inferences  could  be  made  about  the  com¬ 
ponents  of  the  H  vector  in  the  special  case  of  curvilinear  regression.  This  cannot  be  demonstrated 
in  the  general  case.  To  explain  this,  the  pertinent  orthogonal  polynomials  will  be  reviewed. 

In  the  general  case.  Gram-Schmidt  construction  leads  to  the  following  relationships 
(\  l  d  I  1  xi  I 

Pi2  d  1  2  Xi  I  '  ^22  Xi2  <34) 

f\l  'Vt  x  i  I  v  J23  V iS  ‘  dV\  V  ( :i 

and  it  is  seen  that  each  polvnominal  is  a  linear  combination  of  the  independent  variables. 

The  difficulty  lies  in  the  fact  that  generally.  there  is.  no  physical  significance  in  the 
particular  linear  combinations  that  arise.  This  is  more  obvious  when  one  notes  that  these  com¬ 
binations  arise  strictly  according  to  the  manner  in  which  the  independent  variable  values  were 
chosen  and  depend  in  no  way  on  the  functional  relation  between  the  dependent  variable  and  the 
independent  variables. 

Phis  statement  may  be  verified  by  observing  that  the  I)  matrix  is  formed  independently  of 
the  )  matrix. 
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Fortunately,  when  the  Gram-Schmidt  construction  is  applied  to  the  curvilinear  regression 
model,  each  polynominal  is  one  degree  higher  than  the  one  preceding  it 


Pit 

'*11 

Pi  2 

*li 

f  d22  xi 

(37) 

Pis 

'  d\3 

f  d23  Xi  *  rf33  V 

The  nature  of  an  nth  degree  polynomial  is  not  altered  be  its  having  terms  of  degree  less  than  n: 
therefore,  useful  inferences  may  be  made. 

It  is  so  useful  to  form  inferences,  and  to  work  in  the  undistorted  P  space,  that  the  general 
case  will  be  examined  further.  If  only  diagonal  elements  appeared  in  the  l)  matrix,  Eq.  34  could 
be  rewritten 

Pit  xi  i 

P(2  xi  1  *  ^22  X i 2 

Pit  *  (0)  X i  2  *  rf3  3  r,  3 

In  this  speculative  case,  each  p  t  is  related  to  onlv  one  independent  variable.  Therefore,  it 
would  not  be  unreasonable  to  make  inferences  about  the  corresponding  bj  components,  namely, 
that  b.  would  show  a  functional  relation  between  the  dependent  variable  and  the  yth  independent 
variable. 

When  would  this  speculative  case  arise?  The  off-diagonal  terms  in  the  D  matrix  appear 
during  the  Gram-Schmidt  construction.  I'he  only  case  where  they  are  zero  is  the  one  where  the 
columns  of  the  A  matrix  are  already  orthogonal.  This  is  unfortunate,  because  an  orthogonal 
relationship  among  the  independent  variables  is  far  from  a  general  case- 

However,  the  nature  of  an  independent  variable  is  that  its  values  are  free  to  be  chosen 
(hence  the  term  independent).  By  choosing  these  values  so  that  the  A  matrix  is  orthogonal  before 
the  experiment  is  run.  one  gets  a  designed  experiment.  llie  I)  matrix  can  now  be  viewed  as  a 
description  of  the  arbitrary  wav  in  which  values  are  assigned  to  the  independent  variables. 

hor  example,  suppose  that  there  are  two  brands  of  soap,  brand  P  and  brand  ().  Suppose, 
further,  that  each  manufacturer  has  recently  introduced  a  new  version  of  his  brand  containing 
the  ingredient  IBID.  There  are  two  independent  variables:  the  brand,  and  the  presence  or 
absence  of  TRIIJ.  The  value  *■  I  can  be  assigned  to  the  first  independent  variable  (.t  .  ()  when  the 
soup  is  brand  P,  and  the  value  -  I  when  the  soap  is  brand  0  Furthermore,  the  value  +1  can  be 
assigned  to  the  second  independent  variable  (x^j)  when  the  product  contains  THID  and,  -1  when 
it  does  not. 

Io  determine  the  relutive  cleaning  power  of  the  four  soaps,  a  test  agency  meusures  the 
dependent  variable  by  measuring  the  percentage  of  light  reflected  from  sheets  washed  in  the 
various  products.  It  is  desirable  to  compute  statistics,  h from  which  proper  inference  can  be 
made  about  the  two  independent  variables. 


r> 


Hut  not  necessarily  normal 
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Note  the  X  matrix  formed  when  the  following  four  data  points  arc  chosen 

1.  Brand  P  without  TRID 

2.  Brand  Q  without  TRID 

3.  Brand  P  with  TIUI) 

4.  Brand  Q  with  TRID 


The  two  columns  of  ,\  are  already  orthogonal,  and  P  and  D  are 


(51) 


due  entirely  to  normalizing.  Hence 


Components  b  j  P  (  •  }  and  l>  2  P  •  )  may  be  inferred  to  be  die  variation  of  )  with 

brand  name  and  with  TRII).  respectively  This  is  possible  because  a  good  choice  of  data  points 
was  made 

The  p  vectors  in  this  case  are  called  “orthogonal  contrasts."  This  is  only  a  simple  ex¬ 
ample  of  the  designed  experiment  A  specific  class  of  designed  experiments  is  discussed  in  the 
next  section. 

FACTORIAL  ANALYSIS 


One  experimental  design,  leading  to  orthogonal  relationships  among  the  independent  variables, 
is  the  lattice*6  design.  The  analysis  of  this  design  will  be  presented,  and  an  example  will  be  fol¬ 
lowed  throughout  the  section. 

Consider  an  experiment  where  there  are  two  independent  variables,  .t |  and  t2  In  this  example, 
it  is  desirable  that  *|  take  three  distinct  values  and  that  x2  ,n^e  four  distinct  values.  A  special 
situation,  involving  a  total  of  12  points,  occurs  when  the  four  x2  values  are  repeated  for  each  of 
the  three  .tj  values.  In  this  case,  the  paired  values  taken  by  the  two  independent  variables  may 
be  arranged  into  the  following  lattice. 

t,(I).  *2(1)  x,(l).  *3<2>  *,(1).  v2<3)  *,(1),  *2U) 

v,(2).  x2(l)  *,(2),  «2(2)  t,(2).  t2(3)  *,(2).  *2(»  <-,2) 

*,(3).  x2(l)  *,(3).  r2(2)  v ,(1).  x2(3)  *,(3).  *2(« 

\  uid(*|)  tisei!  special  rase  of  (arterial  analx-os 
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where  x,(l),  *j(2),  and  *i(3)  (which  may  be  denoted  Xj(/,),  1  <  j\  <  3)  are  the  three  values  taken 
by  the  first  independent  variable,  and  *2(72),  1  <  j3  <  4,  are  the  four  values  taken  by  the  second 
indspcndent  variable.  In  such  a  design,  j\  is  called  the  “level”  of  Xj,  and  j2  is  called  the  level 
of  x2.  The  variables  xj  and  x2  are  called  “factors.” 

In  this  experiment,  the  dependent  variable  y  may  be  doubly  subscripted  by  the  level  indices 
of  the  independent  variables.  The  experimental  data,  therefore,  may  be  arranged  into  the  same 
lattice 

*11  7.2  713 

721  722  723 

7si  732  733 


7u 

724 

734 


(53) 


To  analyze  such  n  design,  a  regression  model  must  be  set  up.  The  first  step  in  establishing 
a  model  is  to  "partition”  the  independent  variables.  Although  this  term  is  just  now  introduced  in 
factorial  analysis,  it  refers  to  substitutions  already  described  informally  in  the  sections  entitled 
Orthogonal  Polynomials  and  Curvilinear  Regression,  and  the  Design  of  Experiments.  Formally, 
partitioning  is  the  procedure  whereby  the  expression  y  =  /(x)  is  replaced  by  the  expression 
y  *  f(il(x),  i'*(x),  .  .  .  ,  Jtn(  x) .  In  other  words,  the  expression  *y  is  a  function  of  x*  becomes 
"y  is  a  function  of  a  numbe:  f  functions  of  x."  This  procedure  introduces  no  restrictive  assump¬ 
tions. 

A  new  notation  is  used  in  this  formal  definition  for  reasons  that  will  become  apparent  as  the 
model  is  developed.  The  symbol  refers  to  an  arbitrary  function  of  x.  The  superscript  is  an  in¬ 
dex  which  identifies  the  function;  it  is  not  an  exponent.  The  form  of  these  functions  in  the  case 
of  curvilinear  regression  is 

i*  -  1 

X7  *  *  (54) 

x3»*2 

and  so  on,  where,  in  the  conventional  notation  on  the  right  (not  underscored),  a  superscript  is 
indeed  an  exponent.  This  is  a  common  example  of  partitioning,  but  certainly  not  the  only  exam¬ 
ple.  The  arbitrary  functions  may  take  such  diverse  forms  as 

X1  -  1 

X7  =  log  (x)  (55) 

X3  =  tan  (x) 

and  so  on.  The  basic  point  to  remember  is  that  for  a  finite  number  of  x  values,  the  arbitrary  func¬ 
tion  is  only  a  set  of  functional  values  and  may  be  viewed  as  a  vector  in  linear  vector  space. 

Sometimes  only  the  set  of  functional  values  is  given,  and  there  is  no  analytic  form  to  x*.  An 
example  of  such  a  case  may  be  taken  directly  from  the  section  on  the  design  of  experiments: 

X1  *  +  1  if  soap  ia  brand  P  x7  =  +  1  if  TRID  is  contained 

X1--!  if  soap  is  brand  Q  x2  =  -  1  if  TK1D  is  not  contained 

The  latter  case  is  called  a  qualitative  partition,  and  the  former  coses  (such  as  Eq.  54  and  55) 
are  called  quantitative  partitions  . 
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The  regression  model  for  y,  a  function  of  one  variable  with  partitioning,  may  now  be  written 

y  -  v1  +  +  •  •  -+a^n 

and  a  simple  criterion  for  a  linear  model  may  be  expressed 


jL.ma.  (  a  constant)  1  <  i  <  n  (56) 

dx' 

The  lineer  approach  to  lattice  designs  require*  that  each  independent  variable  (factor)  be 

partitioned  into  as  many  arbitrary  functions  as  there  are  functional  values  (levels).  This  is  a 

controversial  point.  The  nonlinear  approach  presented  by  many  texts  does  not  require  partition* 
ing  (in  fact,  partitioning  is  not  possible  without  reverting,  at  least  part  way,  to  a  linear  approach). 
An  objection  to  the  linear  approach  is  that  meaningful  partitions  cannot  always  be  made.  (Take 
three  brands  of  cigarettes  for  an  independent  variable,  as  an  example.)  The  objection  is  not 
valid,  however,  as  meaningless  partitions  can  always  be  performed  and  a  “departitioning"  per¬ 
formed  at  the  end  of  the  analysis.  This  will  be  discussed  further  after  the  development  of  the 
model  has  been  completed  and  the  analysis  has  been  presented. 

In  the  lattice  design  being  followed,  there  are  two  factors,  one  with  three  levels  and  one 
with  four  levels: 

*,(/,)  1  <  /,  <  3 

*2(/2>  1  <  /2  <  4 

The  first  factor  must  be  partitioned  into  three  functions 

A(/j),  Jt2 1(/j),  A(/,)  1  <  ;,  <  3 

and  the  second  must  be  partitioned  into  four 

jt*  Hj2),  xHjJ,  J£3  2(/2),  X  2(/2)  1  <  /'2  <  4 

to  equal  the  number  of  levels  . 

The  partitioning  must  be  performed  so  that  there  is  linear  independence  among  the  arbitrary 
functions  and  hence  in  their  matrix  of  values  .  The  matrix  of  values  for  the  first  factor  is 


*lH\) 

a2‘(l) 

An 

X  ‘(2) 

l'(2) 

l'(2) 

sl(Z) 

m 

i  '(3) 

As) 

If  *j(l)  =  1,  xj  (2)  =  2,  and  *i(3)  =  3,  and  this  first  factor  is  partitioned  according  to  Eq.  54,  the 
matrix  of  values  to  be  checked  for  linear  dependence  is 

n  i  n 


Once  independent  partitioning  is  accomplished  for  each  factor,  the  regression  model  itself 
may  be  considered.  Fora  lattice  design  of  two  factors,  the  general  form  of  this  model  must  be 


rV2  =  /[X*(/’')’ 
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When  the  partition  device  is  applied,  the  model  becomes 

-v/,/j  /£*'  '(/,)•  ±  '(/,)•  ±  2{/V-  i  2(/V-  aV  JL  2< /2) J 

If  the  requirement  is  made  that  r  be  a  linear  function  of  3  |  when  jr2  is  constant,  the  following 
relationship  may  be  derived  from  Kq.  .16 


d\ 


a.  (.*,)  (a  function  onlv  of  x  )  I  ■-  t  .  ^  3 

I 


Furthermore,  if  the  requirement  is  made  that  v  be  a  lineur  function  of  a 2  when  *1  is  held 
constant,  another  relationship  may  he  derived  from  Kq.56 


<9  v 

— ■ —  a  l  x .)  (a  function  onlx  of  x .)  1  £  1,  <  t 

d±2  2 

The  most  general  solution  to  these  differential  relationships  can  be  shown  to  be 


■  V 


l'2 


a,  ,  S  ‘l/jU1  2< *  «|;JL  'l  / ,  >i.2  '< /2>  • 

•  '(/(It  ■’(/'  ,)  •  a22  J?  ’ly  ,)JL  “</,)  * 

3  .  1  3  .  2., 

*  a  "'/V  •  uv2->“  (/VA  (V  • 


a  1  ^  1  / 1^  ^  i 2^ 

°2, 


t  a|4X  ,(/,)x<2(/2) 

.  «24  *'»/,)*  4i/a> 

3  4.. 

.  a,4x  <y,)x  i/2) 


(58) 


provided  that  one  arbitrary  function  of  each  independent  variable  is  merely  a  constant.  The  prac¬ 
tice  of  setting  j.1  1  and  x*  *  'and  for  that  matter  the  first  arbitrary  function  of  any  factor)  to  a 
constant  value  of  unitv  is  so  common  and  so  necessary  for  this  solution  that  it  will  be  assumed 
throughout  the  rest  of  the  text 

This  is.  finullv.  the  regression  model  for  the  3  *  }  lattice  design  being  used  as  an  example. 
That  it  obeys  the  differential  relations  above  may  be  verified  by  substitution. 

As  a  further  demonstration  of  the  linearilv  of  the  model.  Kq.  58  may  be  expanded  into  12 
linear  equations,  in  12  unknowns,  bv  substituting  in  the  12  lattice  points. (The  data(vy|/2  from 
Kq  .53)  are  directly  substituted,  but  the  values  of  the  independent  variables  (  t  |(/|).  *•_,(  )i>  from 
Kq  .52)  must,  of  course,  be  transformed  according  to  the  partitioning  used  t such  as  in  Kq.  54).  *  I 
I'he  u(),j  may.  therefore,  be  determ  ined  by  the  techniques  of  the  previous  sections  However,  a 
more  powerful  technique  is  possible  and  will  be  described  Inter. 


There  are  some  interesting  and  desirable  properties  of  this  model  (Kq.  58)  wh'rh  should  be 
noted.  I’he  terms  in  the  equation  mav  be  arranged  into  the  same  lattice  us  the  design  points 
( Kq.  53).  This  propertv  is  a  kev  to  the  analysis  which  will  be  presented.  Another  property  of 


I '  ><•■•  Kq 


fur  .1  -\mliuli*  of  tin's 
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the  model,  and  one  which  makes  it  realistic,  is  that  when  each  factor  is  partitioned  into  zero- 
order,  linear,  quadratic,  etc.,  functions  (Eq.  54),  the  model  takes  the  form  of  a  Taylor’s  ex¬ 
pansion  of  y  in  two  variables 

r<*,.  *5)  - 


+  **12*2 


+  Q  X 

*  13  *2 


+  0,4*2 


+  02,*,  +oJ2*,*2  +o23  *,*22 


+  °31  V  +  **32  *1*  x2  +  **33  *1*  V  + 


Here  the  level  indices  have  been  stripped  to  enhance  the  form  of  the  expansion,  and  o^  .  is  an 
approximation  for 


MV  VJ  ‘  y <o.6> 

\*xi)  v*2/ 


Because  of  the  finite  number  of  points  in  the  lattice,  only  the  low  order  terms  of  the  expansion 
appear. 

The  model  given  may  be  simplified  (at  least  in  appearance)  as  follows: 

*/,/,  -'*/!  A‘,‘,  w> 

where  the  underscored  capital  letters  are  tensors. 

Tensor  notation  (Ref. 6)  is  intoduced  at  this  point  because  it  is  particularly  suited  to  the 
description  of  the  analysis.  The  small  amount  of  tensor  notation  used  will  be  fully  explained. 

The  four  variables  in  Eq.  60  are  two-dimensional  arrays  similar  to  matrices,  except  that  some  of 
the  indices  appear  as  superscripts.  The  elements  of  the  Xj  j  array  appear  in  Eq.  S3. 

The  elements  of  X1  are*  *(/.),  where18 
'  t 

1 1  <  <  3 

and  those  of  X  :2  are.£  2(/' ),  where 
'2 

1  <  *2  <  4 
1  <  i2  <  4 

The  elements  of  the  X  ;  tensor  span  the  a  coefficient  of  Eq .  58. 

1*2 

Tensors  may  be  multiplied  in  two  different  ways  to  obtain  either  an  inner  or  an  outer  product. 
Only  the  inner  product  is  of  interest  in  this  development,  and  only  that  product  will  appear. 

As  a  first  example  of  the  tensor  product,  consider  the  product  X-1  X2  from  Eq.  60.  This 

/  I  1  2 
l  (  *  C 

product  yields  a  four-dimensional  array  which  may  be  denoted  X1 .2  without  notations!  ambiguity. 

/j/2 

18 

"These  elements  may  be  recognized  as  the  matrix  of  values  (Eq.  57). 
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The  /jij/j/j  th  element  of  this  array  is  the  product  of  the  ij/j  th  element  from  Xj^  and  the  »2/2 th 
element  from  Jy  1  . 

The  continuation  of  the  product  in  Eq.  60  presents  a  second  example  of  the  tensor  product 

±. . 

i ih  *1*2 


Note  that  the  superscripts  of  X  appear  as  the  subscripts  of  X- 

There  are  two  steps  in  forming  this  product.  The  first  is  forming  the  following  four¬ 
dimensional  array  (there  are  six  indices,  but  only  four  are  distinct): 

yV2. 

~ V2V2 

where  the  *\/[,272  th  element  of  this  array  is  formed  from  the  product  of  the  *|/j*2/2th  element 
of  X  and  the  »ji,  th  element  of_j4. 

The  inner  product  is  completed  by  a  second  step  where  the  indices  appearing  both  as  subscripts 
and  as  supeiacripts  are  to  be  summed  out  as  follows: 


‘/l'2 


=  1  X  v'l'a.  . 
‘l  »2  V2*l‘2 


(61) 


which  reduces  the  product  to  a  two-dimensional  array.  A  numerical  example  of  such  a  tensor 
product  appears  in  Appendix  E. 

Returning  to  the  subject  of  factorial  analysis,  in  Eq.  60  there  appears  a  tensor  for  each  in¬ 
dependent  variable  or  factor ■  JC . 1  and  X,2  •  While  their  product  can  be  immediately  written 

/ 1  '2 
i  i 

X:‘  2 ,  it  is  better  to  analyze  them  separately.  Each  X  tensor  can  be  analyzed  in  a  manner 
•  li  2 


as 


similar  to  the  matrix  methods  given  previously. 

It  is  desirable  to  separately  transform  each  of  the  X  tensors  into  orthonormal  £  tensors. 
This  may  be  done  by  the  Gram-Schmidt  construction,  which  was  previously  described  by  the 
matrix  transformation 

P=  X  D 


This  transformation  may  be  rewritten  in  tensor  notation 

£sj  =  Xj  11*  (62) 

Some  insight  into  tensors  can  be  gained  by  studying  this  equation.  Because  it  describes  the 
same  transformation  as  the  matrix  equation,  it  -  an  be  observed  that  the  superscript  of  the  ten¬ 
sor  embodies  the  row  property  of  the  matrix*®  and  the  subscript  of  the  .ensor  embodies  the 
column  property.20 

Consider  now  a  tensor  Cj,  such  that  the  tensor  product 

£?£*=4 

where  the  elements  of  J.  are  unity  if  i  =  k  and  are  zero  if  i  4  k.  Such  a  tensor  £,  contains  the  in¬ 
verse  elements  of  the  H  array  in  the  matrix  sense. 

Also  called  the  contravariant  property. 

20  Also  called  the  covariant  property. 
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The  tensor  i  has  the  property  that,  when  it  multiplies  another  tensor  such  asX  the  follow¬ 
ing  relation  results: 

Performing  this  multiplication  according  to  the  inner  product  rules  will  show  that  the  individual 
elements  of  X1-  ore  equal  to  the  individual  elements  of  Jl|,  so  the  change  is  only  symbolic. 

This  relationship  will  be  used  below. 

Both  sides  of  Eq.  62  may  be  multiplied  by  £  to  give  the  following: 

ISjLk.  =  2!-£?£*  =■!•  (63) 

A  further  multiplication  by  £k  gives 

(64) 

This  relation,  with  the  indices  properly  subscripted,  may  be  substituted  into  Eq.  60  for  each 
factor  X  to  give 

I  ■  =£'ldl  E'1  £<2  4  •  (65) 

h>2  1 1  N  >2  *2  ‘l‘2 

Some  desirable  properties  of  tensors  allow  a  rapid  completion  of  the  analysis.  Tensors,  un¬ 
like  matrices,  obey  the  commutative  rule  as  well  as  the  associative  rule  of  multiplication;  for 
example 

f ;  Cj  =  E]  (66) 

Therefore,  one  may  rearrange  the  product  terms  in  Eq.  65  at  will.  One  way  follows: 


where 


Y  •  =/■'  /2  B 

~V  2  '2~V 


It  is  epeated  at  this  point  that  each  £.  tensor  is  a  two-dimensional  array  which  contains  the 
same  elements  as  an  orthonormal  matrix  P.  \  property  of  such  matrices  is  that  Pf  P  -  /.  This 
relation  may  also  be  rewritten  in  tensor  notation 

£[  ESj  =lat 

It  is  legitimate  here  to  use  the  same  s/mbol  for  each_£  tensor,  because  each  contains  the  same 
elements.  However,  the  index  which  describes  the  column  property^*  in  the  original  Ps  array 
must  describe  the  row  property  in  the  new  P i  array.  1 

Hence  j  is  a  subscript  in  the  tensor  corresponding  to  the  matrix  P,  and  is  a  superscript  in  the 
tensor  corresponding  to  the  matrix  P /  It  is  in  this  way  that  tensors  are  transposed.  Both  sides 


2 1 

1  This  is  because  the  matrix  P  was  defined  so  that  it  was  column-orthogonal. 
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of  Eq.  67  may  be  multiplied  by  and^2  to  give  the  following:22 


D*1  n*  2  y  _  d 
X‘l  ‘  2  * / |/' 2  "  "‘1*3 


(68) 


The  tensor  4  may  be  solved  for  in  a  manner  analogous  to  the  matrix  solution  of  Eq.  25. 
However,  as  before,  this  will  be  considered  an  anticlimax  in  the  analysis,  and,  in  fact,  it  is  very 
seldom  done  in  factorial  analyis. 

Equation  68  may  be  compressed  by  the  rule  of  tensor  products  as  follows: 


>Va 


“*  i*2  i \i 2 


a,  , 

*  1*2 


(69) 


The  tensor  operates  on  the  data  tensor  2’, to  yield  a  third  tensor  B.  A  transformation  is 
presented  in  the  next  section  which  changes  Jf  and  B  into  one-dimensional  arrays  while  pre¬ 
serving  the  total  number  of  elements  contained.  The  same  transformation,  when  applied  to  the 
P  tensor,  reduces  it  to  a  two-dimensional  array.  It  will  be  shown,  further,  that  this  latter  array 
is  also  orthonormal.  This  will  establish  that  the  elements  of  the  ft  tensor  are  components  of  the 
lata  2’.  Therefore,  the  theorem  of  Pythagoras  applies,  which,  for  the  example,  is  expressed 


where  a  double  summation  is  required  because  of  the  double  index. 


(70) 


Furthermore,  the  elements  of  B  may  be  inferred  to  be  variations  of  the  data  due  to  certain 
/ariations  in  the  independent  variables.  To  complete  the  analysis,  it  is  necessary  to  identify 
hese  variations.  The  key  to  this  identification  is  the  double  subscript  of  the  ft  tensor.  Just 
is  the  double  subscript  of  the  tensor  identifies  the  individual  data  with  the  levels  taken  by 
he  independent  variables,  the  subscripts  of  the  B  tensor  (which  run  the  same  gamut  as  the  Y 
lubscripts)  identify  the  components  with  these  variations  in  the  independent  variables. 

A  study  of  this  can  be  made  by  expanding  Eq.  67  according  to  the  definition  of  the  tensor 
>roduct 


4iiP  ‘(/V  P  2<h)  +  bl2  P 
k 21  P  l(/l}£  2}  +  -  2  2  jP  ^ 

3  .  1  -  3 

-  3  1  P  (/V  P  W  +  -32  P 


1(/V £J2(/2)  +£i3p,l(/Vp  20'2)  +  £i«P  1(/Vp  2</2)  + 

‘t/V  +  ^23f*I(/l)  f3^/  2)  +  t2ip2^ii)  + 

*(/,  )£**(/ 2)  +  ^33f  3  1(/  l>  P  3  *</ 2>  +  ^34f3l(/1)p<2</2) 


(71) 


22 

The  difference  between 
umbers.  12 


m6  2 


is  only  symbolic;  they  each  contain  the  same  array  of 
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As  an  example,  let  the  independent  variables  be  partitioned  according  to  Eq.  54. 


2* *(/ 1>  “  1  A/V"1 

*,‘(/l)-x1  s2(jj-x  2 

£^i2)ss22 
x*^i2)  =  x23 

Then  from  Eq.  62 

f'v-rfl; 

£  'Vil-rfjJ  +  rfa,1*! 


p  2(/2)=i,2 
*“  2 

p\^i\]+4\\x2 

£%>  =  A 'A  +  A’  +  A* 


11  2  2 
Therefore,  p  1  and  p  2  are  zero-order  functions;  p  1  and  p  2  are  linear  functions  of  their 

3  3  4  2 

respective  factors;  p  1  and  p  2  are  quadratic  functions  of  their  respective  factors;  and  p  is  a 

cubic  function  of  x2- 

An  inspection  of  Eq.  71  will  show  that  an  element  b ,  of  the  B  tensor  shows  the  variation 

'sls2 

in  the  data  Y,  due  to  the  variation  in  p  1  and  in  p  2.  A  component  which  depends  on  two  varia¬ 
tions  is  called  an  interaction  component.  In  a  completely  general  case,  all  bs  s  are  interaction 

components.  However,  in  practical  cases,  such  as  the  example  being  given,  the  independent 

variables  are  partitioned  in  such  a  way  that  the  first  orthogonal  function,  p1  of  each  is  zero  order, 
or  constant.  Because  of  this,  certain  degenerations  occur.  The  component  bll  is  inferred  to  be 
the  variation  of  the  data  due  to  the  zero-order  variation  of  Xj  and  of  x  2-  This  interaction  is  doubly 
degenerate  and  is  itself  zero  order.  The  term  £,2  shows  the  variation  of  the  data  resulting  from 
the  zero-order  variation  of  *]  and  the  linear  variation  of  *  2>  This  interaction  is  singly  degenerate 
and  is  called  a  main  effect,  namely,  the  linear  x  2  effect.  The  rest  of  the  inferences  are  found  in 
the  same  manner  and  are  tabulated  in  Table  1. 

Henceforth,  only  interactions  which  do  not  degenerate  will  be  termed  interactions. 

An  inference  need  not  always  be  made  about  a  component.  In  order  to  test  whether  components 
are  showing  a  true  functional  relation  between  y  and  the  independent  variables,  it  is  necessary  to 
have  extra  components  where  no  functional  relation  is  expected  to  appear.  These  components  are 
of  the  same  nature  as  those  described  in  the  section  on  curvilinear  regression  and  are  termed  error 
components.  An  F  test  can  be  made  where  each  component,  about  which  an  inference  is  made,  is 
tested  against  the  error  components. 

In  the  example  given,  there  are  not  enough  components  to  test  all  possible  inferences.  Two 
approaches  are  possible,  hirst,  it  may  be  that  there  are  no  true  functional  relations  expected 
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TABLE  1.  Inferences  for  a  3  x  4  Experiment  With  Linear,  Quadratic, 


and  Cubic  Partitioning 


Component 

x  |  Variation 

x2  Variation 

Inferences 

Sit 

zero  order 

zero  order 

zero  order0 

fi  12 

zero  order 

linear 

lin.  x2 effect 

« 13 

zero  order 

quadratic 

quad.  *2  effect* 

fiu 

zero  order 

cubic 

cub.  x2  effect4 

Sjl 

linear 

zero  order 

lin.  x  1  effect4 

fi  22 

linear 

linear 

lin.  xj  lin.  x2  interaction 

fi  23 

linear 

quadratic 

lin.xi  quad.  x2  interaction 

fi24 

linear 

cubic 

lin.  x  j  cub.  x2  interaction 

B3I 

quadratic 

zero  order 

quad,  x  j  effect4 

B32 

quadratic 

linear 

quad.  Xjlin.  x2  interaction 

fl33 

quadratic 

quadratic 

quad.  X]  quad.  x2  interaction 

«  34 

quadratic 

cubic 

quad.  Xj  cub.  x2  interaction 

“Doubly  degenerate, 
4  Singly  degenerate. 


among  the  interaction  components.  In  this  case,  the  F  ratio  used  to  test  a  true  functional  relation 
between  y  and  a  linear  variation  in  Xj  is  computed  as  follows:23 

F(l,6)  — -  (*21> 


+^23*  +  (^24)  +  ^32^  +(^33)  +  ($34)] 

and  for  a  linear  relation  between  y  and  x 2 


(72) 


F(l,6)  = 


<412> 


2  2  2  2  2  2 

6„„)  +^23^  +  +  (k,i)  +  (£.,4) 


2  2' 


-2*' 


:  32' 


-33' 


•3* 


(73) 


and  so  on  for  quadratic  and  cubic  relations.  If  it  is  necessary  to  infer  that  interaction  effects  do 
exist,  then  it  is  necessary  to  create  more  components,  which  is  the  second  approach.  The  only 
way  to  do  this  is  to  add  a  new  factor.  This  factor  may  be  some  new  variable  in  the  analysis,  or 
it  may  be  a  factor  about  which  no  inference  is  made.  The  latter  case  is  called  multiple  replication, 
because  it  is  merely  repetition  of  the  experiment. 

It  is  not  always  possible  to  partition  a  factor  in  a  meaningful  manner.  A  “brand  name”  factor 
was  given  as  an  example.  If  factor  x  j  of  that  example  had  represented  brands  A,  B,  and  C  of 
automobile  tires,  then  it  would  be  meaningless  to  talk  about  the  linear  and  quadratic  variation 
among  the  brands.  However,  one  can  still  make  the  meaningless  partition  and  get  the  total  varia¬ 
tion  of  y,  due  to  brands  by  adding  (vectorially)  the  two  meaningless  components 


—  21  *  —  3  1  =  ^brands,  1 


Because  b21  and  b3l  are  components  in  orthogonal  directions,  the  length  of  the  new  vector  is 


(b2l)2  *  (b 3 j)2 


by  the  Pythagorean  theorem.  This  may  be  called  departitioning. 

23 

Note  that  the  denominator  contains  the  interaction  components  shown  in  Table  1. 
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If  one  uiumei,  ae  before,  that  the  interaction  terma  are  error,  the  F  ratio  to  teat  the  varia- 
tion  of  y  due  to  branda  ia 


F(  2,6) 


_ T  [<»«>* +<4n>a] _ 

+^33>  +(ia4)‘ +  (A,3>  +  (i33)  +  <£,,)] 


where  the  general  rule  ia  applied  that  the  F  function  ia  the  ratio  of  the  mean  square  components 
being  tested  to  the  meon  square  components  considered  to  be  due  to  experimental  error. 

The  choosing  of  the  components  to  be  tested,  and  of  the  F  ratios  to  be  formed,  is  called  the 
summary  of  the  analysis  and,  as  stated  before,  is  strictly  within  the  working  area  of  the  statisti¬ 
cian. 

Thus  far,  only  a  particular  two-factor  experiment  has  been  considered.  Generalizing  to  the 
n  factor  experiment  ia  remarkably  simple  in  tensor  notation  and  is  summarized  below. 

SUMMARY  OF  THE  ANALYSIS  OF  LATTICE  DESIGNS 

Assume  the  regression  model 

v  _  1  v‘3  i  (75) 


1.  Partition  each  independent  variable  to  obtain  the  matrix  of  values  which  are  to  be  inserted 
into  the  individual  X  tensors. 


Sl(j  j)  *  j/H*  ,(/',)! 
*<2(/2)  =  2(/  2>] 
x  *(/3)  =  & 


(each  x  1  an  arbitrary  function  of  *j) 
(each  J  2  an  arbitrary  function  of  x2) 
(each  s  3  an  arbitrary  function  of  x3) 


where,  for  each  factor,  the  number  of  values  taken  by  i  (number  of  functions)  is  equal  to  the 
number  of  values  taken  by  j  (number  of  levels)  to  obtain  a  square  matrix  of  values. 

2.  If  each  matrix  of  values  possesses  linear  independence,  form  the  orthonormal  tensors  by 
the  Gram-Schmidt  construction 

1 1  *  1  *  1 

jP'2  =  j‘2d'2 

/  2  >2  2 

P*3  =  A'3  D*3 

I  3  >  3  3 

3.  Transpose  each  £  and,  using  the  solution  process  of  Eq.  67  and  Eq.  68,  evaluate  the 
following  equation 


-*152,3 


p*l*2*3’  ’  •  y 

■  ’  1  •-/,/ j/,*  •  * 
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4.  Form  inferences  about  §  (which  may  now  contain  three-way  and  higher  interactions)  and 
summarize  according  to  the  desires  of  the  statistician. 

The  factorial  analysis  program  in  Appendix  C  uses  this  method.  Other  facets  of  this  program 
will  appear  in  succeeding  sections. 

The  reader  is  referred  to  Ref.  7  for  another  description  of  the  linear  approach  to  lattice 
designs. 


TRANSFORMATION  OF  INDICES 

As  stated  in  the  preceding  section,  a  transformation  will  be  described  which  reduces  the 
multidimensional  Y  and  B  arrays  to  one-dimensional  arrays.  The  transformation  works  on  the 
indices  of  the  tensor,  but  it  does  not  alter  the  number  of  elements  in  the  array  or  their  numerical 
values. 

This  transformation  is  of  more  than  academic  interest,  because  the  addressable  memory  of 
a  digital  computer  is  one-dimensional.  So,  an  n-dimensional  array  must  be  represented  by  a  one¬ 
dimensional  array,  the  elements  of  which  are  found  by  a  single  subscript  computed  from  the  many 
subscripts  of  the  n-dimensional  array.  ' 

Assume,  for  example,  that  there  are  m  factors  with  the  number  of  levels  denoted  by  f  j,  f2, 

....  /  .  The  total  number  of  data  points,  n,  is  equal  to  the  product  of  the  fa.  Consider  the 

data  tensor  Y,  .....  ,  a  variable  with  m  subscripts.  A  single  subscript  j  may  be  computed  from 
“'1)2  'm 

these  subscripts  as  follows:^*  1 

i- 1  =  /,-! +  +  (/3-iVi  f2  +  </4-iV,  f2  f3  *  ■  •  •  +  (/,-D/i  f 2 '  ‘ '  fm—  i  <77> 

The  above  formula^  will  be  recognized  immediately  by  one  who  has  coded  operations  on 
m-dimensional  arrays  in  digital  computing  machine  language,  and  the  number  theorist  will 
recognize  it  as  a  representation  of  /-I  in  a  hybrid-base  number  system.  That  is,  a  number  may  be 
written  (/m-l)  •  •  •  ( y 2—  1 )  (/'  j—  1 ),  which  does  not  represent  /-l  in  a  decimal  or  a  binary  number 

system,  but  in  a  system  where  the  low  order  digit,  (/j-1),  is  of  base  /jjlhe  second  order  digit,  ( / 2—  1 ) , 
is  of  base  f 2,  and  so  on.  Because  the  last  viewpoint  delivers  the  greatest  insight  into  the  trans¬ 
formation,  an  example  will  be  given. 

Consider  a  three-factor  experiment  where  factor  Xj  has  two  levels,  *2  has  three  levels,  and 
x3  has  four  levels.  Table  2  illustrates  how  the  three-dimensional  array  may  be  “stretched  out* 
into  a  single  array.  An  important  property  of  the  transformation  in  Table  2  is  that  it  is  a  one-to-one 
transformation.  That  is,  for  every  j  there  is  a  unique  set,  j v  jT  ,  jm,  and  vice  versa.  The 
uniqueness  may  be  proved  by  induction  in  a  manner  similar  to  the  proof  in  Ref.  8. 

^This  j  is  notationally  distinct  because  it  is  not  subscripted. 

^A  cumbersome  translation,  j  —/-J,  is  used  because  the  FORTRAN  language  does  not  allow  indices  to 
run  from  zero. 
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TABLE  2.  Transformation  of  Indices  fora 
2x3x4  Experiment  Showing  the  Base — 10 


and  Hybrid— Base  Representation  of  the 
Single  Index  Formed 


Jr\ 

JT\ 

/-I 

(Base  10) 

M 

(Hybrid  base) 

0 

0 

0 

0 

000 

1 

0 

0 

1 

001 

0 

1 

0 

2 

010 

1 

1 

0 

3 

Oil 

0 

2 

0 

4 

020 

1 

2 

0 

5 

021 

0 

0 

1 

6 

100 

1 

0 

1 

7 

101 

0 

1 

1 

8 

110 

1 

1 

I 

9 

111 

0 

2 

1 

10 

120 

1 

2 

1 

11 

121 

0 

0 

2 

12 

200 

1 

0 

2 

13 

201 

0 

1 

2 

14 

210 

t 

1 

2 

15 

211 

0 

2 

2 

16 

220 

1 

2 

2 

17 

221 

0 

0 

3 

18 

300 

1 

0 

3 

19 

301 

0 

1 

3 

20 

310 

1 

1 

3 

21 

311 

0 

2 

3 

22 

320 

1 

2 

3 

23 

321  ’ 

This  transformation  may  be  applied  to  the  subscripts  j  and  t  in  Eq.  69  as  follows: 

flXj-£t  (78) 

It  is  important  to  recognize  that  Eq.  78  describes  the  same  arithmetic  computations  as  Eq.  69. 

In  fact,  the  computer  programmer  will  recognize  Eq.  78  as  the  method  of  coding  the  computations 
of  Eq.  69  in  machine  language. 

Equation  78  may  be  immediately  rewritten  in  matrix  notation 

P'Y  =  B  (79) 

noting  that  the  simplified  tensor  products  of  Eq.  78  corresponds  exactly  to  a  matrix  product. 
Equation  67  may  be  transformed  in  the  same  manner 

y  =  P  B  (80) 

applying  first  the  rule  of  tensor  products  to  obtain  Y  ■  ■  =  P*1*  2  B 

~>1>2  ~lll2~sls2 
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The  only  way  in  which  Fq.  79  and  80  may  both  be  true  is  that 

P'P  =  I  (81) 

Thus  it  is  shown  that  die  matrix  P.  formed  by  the  transformation  of  subscripts  and  superscripts 


-  ,  n 4  1 5  2  ■  .  , 

of  the  tensor  P ■  .  ,  is  orthonormal  It  is  therefore  possible  to  observe  two  things.  First,  the 
/  2 

lattice  design  fulfills  Jie  desired  requirement  that  the  independent  variables  be  mutually 
independent  (orthogonal).  Second,  the  bs  s  _  are  truly  components  of  the  data^-  y  .  and 

the  inferences  of  the  previous  section  are  correct. 


FRACTIONAL  REPLICATION 

The  lattice  designs  already  described  are  not  the  only  means  by  which  an  orthogonal  rela¬ 
tionship  may  be  established  among  the  independent  variables  of  an  experiment.  It  is  not  within 
the  scope  of  this  report  to  describe  the  others  in  great  detail.  One  type,  however,  is  so  common, 
and  so  intimately  related  to  lattice  designs,  that  something  should  be  said  about  what  it  is,  why 
it  is  used,  and  how  it  may  be  handled  on  digital  computers. 

Fractional  replication  is  related  to  what  is  called  the  redundancy  check.  The  reader  will 
probably  be  familiar  with  the  parity  check  on  binary  numbers  or  with  check  sums  on  punched 
cards.  The  redundancy  check  may  be  described  as  taking  essential  information  and  including 
with  it  redundant  information  computed  from  the  essential  information  in  some  predetermined  way. 
The  purpose  of  this,  of  course,  >s  to  detect  transmission  errors  (or,  when  the  human  element  is 
involved,  to  reduce  the  possibility  of  forgery  of  numerical  information).  One  may,  for  instance, 
look  at  the  10-digit  serial  number  on  his  oil  company  card  and  note  that  the  company  could  not 
possibly  have  10  billion  credit  customers  as  this  number  implies.  If  some  of  these  numbers  are 
used  as  a  redundancy  check,  one's  chances  of  forging  a  legitimate  serial  number  are  small  .  A 
redundancy  check  may  be  added  to  the  subscripts  of  data  tensors  from  lattice  designs.  A  par¬ 
ticular  class  of  redundancy  checks  will  be  described  which  requires  that  two  or  more  factors 
have  the  same  number  of  levels.  For  example,  consider  a  simple  two-factor  experiment  where 
each  factor  has  three  levels.  The  Y  tensor  than  has  nine  elements 

^11.^12*  ^13’  y 2  I  *  ^ 2 2’  ?2  3’  ^31’  ^32’  ?33 

A  third  subscript  may  be  added  which  has  the  property  that  the  sum  of  the  subscripts  modulo 
3  is  zero 


yHI’^123’  ^  l  32'  ^21  3’  ?22  2’  ^231’  ?3!2’  ^321'  ?3  3  3 

Consider  now  the  case  where  a  new  three-level  variable  is  introduced  which  takes  the  level  sug¬ 
gested  by  the  third  subscript.  In  this  case,  some  interesting  things  happen.  First,  the  third 
subscript  is  now  a  redundancy  check  in  mathematical  form  only.  Actually,  it  achieves  the  same 
status  as  the  other  two  subscripts.  Second,  a  three-factor  experiment  has  resulted,  but  some¬ 
thing  is  missing,  namely,  those  elements  of  a  27-element  data  tensor,  the  sums  of  whose  sub¬ 
scripts  equal  1  and  2  modulo  3.  In  fact,  there  is  one-third  of  a  3  x  3  x  3  experiment  (written 
1/3  x  33).  This  is  called  a  fractional  replicate  of  the  full  33  experiment. 
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The  three  subscripts  possess  the  property  of  symmetry.  'Phut  is,  if  originally  there  were  only 
the  second  and  third  subscripts,  the  value  of  the  first  would  be  determined  bv  the  same  rule.  In 
fact,  if  the  first  subscript  is  dropped,  the  second  and  third  take  the  same  set  of  vulues  that  the 
first  and  second  did  originally  (but  not  in  the  same  order);  that  is,  all  the  points  on  a  3  x  3 
lattice  appear.  Similarly,  the  second  subscript  is  a  redundancy  check  on  the  first  and  third,  and 
if  this  second  subscript  is  dropped,  the  first  and  third  contain  all  the  points  on  a  3  x  3  luttice 
(and  again,  the  order  is  different). 

Becuuse  it  is  a  lattice  design,  the  original  two  factors  have  an  orthogonal  relationship  be¬ 
tween  them.  But,  bv  the  property  of  symmetry,  it  must  be  concluded  that  the  third  is  orthogonal 
to  each  of  the  original  two.  Thus,  the  fractional  replicate  is  a  designed  experiment. 


However,  a  paradox  exists.  Although  there  are  three  subscripts,  there  is  onlv  a  two* 
dimensii  nul  urray .  In  fact,  the  analysis  made  for  complete  factorial  designs  may  be  made  here, 
ignoring  any  one  of  the  three  subscripts  (because,  as  already  noted,  in  all  three  cases  the  re¬ 
maining  subscripts  describe  a  complete  lattice).  But  each  process  analyzes  only  two  factors. 
What  happens  to  the  third/ 


The  solution  tensor  U  consists  of  nine  elements;  the  zero-order  components  of  the  data,  the 
two  main  effects  of  the  first  factor,  the  two  main  effects  of  the  second  factor,  and  four  interaction 
components.  Because  die  vector  basis  is  complete,  the  components  related  to  the  third  fuctor 
must  be  contained  somewhere.  Bv  its  very  nature,  the  zero-order  component  is  eliminated,  and, 
as  alreadv  noted,  the  thiee  main  factors  are  niutuallv  orthogonal.  Therefore,  the  two  components 
of  the  third  factor  mast  lie  along  some  of  the  interaction  components  of  the  other  two  factors. 


In  this  case,  each  factor  is  said  to  be  "confounded”  with  the  interaction  of  the  other  two 
factors  This  dual  identity  of  c  omponents  presents  a  problem.  How  does  one  know  whether  a 
component  results  from  one  factor  or  the  interaction  of  two  others/  Proper  inference  can  solve 
the  pioblem.One  muv  have  good  reason  to  believe  that  there  are  really  no  interaction  effects. 
He  then  infers  the  components  to  be  onlv  main  effects 


One  wav  to  have  confidence  in  such  inferences  is  to  confound  fac  tors  with  high-order  inter¬ 
actions.  (In  the  previous  section,  the  regression  model  was  described  as  an  estimate  of  luylor  s 
expansion  of  the  dependent  variable.  One  therefore  expects  components  resulting  from  high-order 
effects,  even  when  real,  to  be  small,  because  of  the  reciprocal  factorials  in  the  high-order  terms 
of  the  expansion.)  The  onlv  way  to  confound  high-order  interactions,  in  the  class  of  fractional 
replication  being  described,  is  to  have  many  factors  with  the  same  number  of  leve  ls.  Consider,  for 
instance,  a  I  x  t  ■  I  ■  I  experiment.  The  independent  variable  will  have  four  subscripts,  each 
taking  values  I  to  (.  One  may  add  a  fifth  subscript  with  a  redundance  check,  so  that  the  sum 
of  the  five-  must  be  0  modulo  I.  and  may  thus  introduce  a  fifth  factor,  f.ach  factor  will  be  con¬ 
founded  with  a  four-way  interaction.  One  then  infers  confidently  that  the  components  show 
onlv  the  functional  relation  of  the  factor  and  not  a  four-way  interaction. 

This  example  also  illustrates  the  principal  advantage  of  fractional  replication.  Kith  four 
factors,  one  alreadv  had  a  need  for  236  experimental  points.  Without  fractionating,  a  fifth  factor 
would  require  1021:  but  with  it,  the"  need  is  still  for  only  236. 
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The  type  of  fractional  replication  described  above  ia  called  the  Latin-square  class,  because 
the  simplest  cases  lead  to  Latin-square  patterns.  For  instance,  from  the  1/3  x  33  example,  one 
may  set  up  the  following  square: 


1 

2 

3 

1 

1 

3 

2 

2 

3 

2 

1 

3 

2 

1 

3 

where  the  first  subscript  appears  across  the  top,  the  second  along  the  side,  and  the  third  in  the 
boxes. 

There  are.  of  course,  more  complex  examples  in  the  Latin-square  class;  for  instance,  a 
1A3  x  4  x  33  experiment.  The  fo  ur-level  factor  has  no  redundancy  check  and  hence  no  con¬ 
founding  with  other  factors.  The  other  three,  considered  separately,  form  the  classical  Latin 
square  above. 

There  are  other  types  of  redundancy  checks  which  give  rise  to  other  factorial  models.  In 
some  cases,  there  can  be  two  independent  redundancy  checks  on  a  number  of  factors  of  the  same 
number  of  levels.  This  is  called  the  Greco-Latin-square  class  of  fractional  replication.  Other 
types  of  checks  exist  which  give  rise  to  incomplete-block  and  other  exotic  designs,  some  of 
which  can  be  solved  by  the  programming  in  appendixes  A,  B,  and  C  and  all  of  which  can  be 
solved  by  the  methods  of  this  report. 

There  are  two  ways  of  analyzing  fractional  replicates  by  the  methods  described  under  Fac¬ 
torial  Analysis.  One  may  be  called  “putting  the  burden  on  the  statistician'  and  the  other  “putting 
the  burden  on  the  computing  machine.'  In  the  first  method,  the  redundant  subscripts  are  ignored, 
and  the  regular  factorial  analysis  is  made,  which  computes  each  single  degree  of  freedom.  The 
statistician,  knowing  the  redundancy  check,  then  searches  among  the  high-order  components  for 
the  missing  factors.  The  second  method  requires  that  the  rest  of  the  lattice,  of  which  there  is  only 
a  fraction,  be  filled  in  with  zeros.  The  regular  factorial  analysis  is  then  made,  and  highest- 
order  interactions  are  successively  discarded  until  there  are  fewer  components  than  data.  The 
remaining  components  are  divided  by  the  fraction  of  the  experiment.  In  this  method,  all  com¬ 
ponents  are  properly  identified  as  described  under  Factorial  Analysis. 

As  an  example  of  the  second  method,  consider  the  1/3  x  33  experiment.  This  is  first  con¬ 
verted  into  a  33  by  filling  those  elements  of  a  3  x  3  x  3  lattice  that  are  not  in  the  design,  with 
zeros.  The  factorial  analysis  is  then  made.  The  two-  and  three-way  interaction  components  are 
discarded,  and  the  main  effects  are  divided  by  1  3.  Care  must  be  taken  in  calculating  F  ratios 
in  this  case.  The  denominator  must  be  calculated  in  the  manner  of  Eq.  46,  where  the  squared 
components  to  be  tested,  including  the  zero-order  component,  are  subtracted  from  the  total  sum 
of  squares  of  the  data. 


THE  MISSING-VALUE  PROBLEM 

It  sometimes  happens  that  one  or  more  points  of  a  desired  experiment  are  missing.  This  can 
happen  for  any  number  of  reasons,  such  as  a  failure  in  recording  instruments  or  an  operator  er¬ 
ror.  What  is  worse,  in  some  cases  these  points  are  irreplaceable. 

When  there  are  points  missing,  there  are  two  possible  approaches  to  analyzing  the  data.  One 
approach  is  to  compute  estimates  of  the  missing  value  on  the  basis  of  the  components  of  the  data 
about  which  real  functional  relationships  are  inferred.  The  other  approach  is  to  revert  to  the  gen¬ 
eral  regression  model. 
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There  are  disadvantages  in  both  approaches,  and  both  are  subject  to  much  criticism.  In  the 
first  method,  there  is  the  disadvantage  that  the  solution  to  the  regression  model  that  is  formed 
after  the  computation  of  missing  values  is  no  longer  the  best  solution  in  the  sense  of  least 
squares.  It  has  the  advantage,  however,  that  the  inferences  are  preserved. 

In  reverting  to  the  general  regression  model,  one  has  the  advantage  of  an  exact  solution  in 
the  sense  of  least  squares,  but  the  inferences  become  distorted  and  F  tests  ore  not  properly 
made  Another  advantage  here  is  that  it  is  a  noniterative  approuch,  whereas  the  missing-vulue 
approach  rcmircs  iteration  if  there  are  two  or  more  missing  values.  On  a  digital  computer,  this 
last  is  not  an  important  consideration 

The  method  of  \illarsUtef.  9)  has  served  the  authors  well  us  a  missing-value  technique  and 
has  been  coded  in  machine  language  for  the  lti\f  70  f  (Hcf.  10).  I’nfortunalely .  it  has  not  been 
updated  to  I'DHTH  VV 


FINAL  REMARKS 

KOtlTHAVcodcd  programs  for  general  regression,  curvilinear  regression,  and  factorial  anal¬ 
ysis.  using  the  methods  of  the  text,  appear  in  Appendixes  A.  M.  and  Although  these  programs 
have  output  summaries,  thev  cannot  be  described  as  final  Tliev  are  presented  muinlv  to  illustrate 
the  I  inear  method.  The  potential  user  is  invited  to  take  them  and  choose  his  own  methods  of  sum¬ 
marizing  the  output  and  of  handling  fractional  replication  and  missing  values. 
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Appendix  A 

USE  OP  THE  MULTILINEAR  REGRESSION  PROGRAM 

The  data  input  may  be  described  best  by  referring  to  the  basic  matrix  equation 

XA  =  Y 

where  X  is  an  n  by  m  matrix,  A  is  a  column  vector  of  m  elements,  and  Y  is  a  column  vector  of  n 
elements.  The  unknown  is  the  A  vector. 

Card  1,  FORMAT  (215,  F5.0,  54H  .  .  .),  has  the  following  parameters: 

N  (1)  =  m  which  is  less  than  51 . 

N  (2)  *  n  which  is  less  than  1001 . 

SF  =  scale  factor  for  F  vector  (power  of  10). 

ID  =  identification  of  problem  with  any  legal  Hollerith  characters, 
m  x  n  -  product  less  than  16,001  . 

Card  2,  FORMAT  (12A6),  contains  identification  for  each  of  the  m  columns  in  the  X  matrix. 
Use  as  many  cards  as  needed. 

Card  3,  FORMAT  (12F6.0),  contains  the  values  of  the  elements  in  the  Y  vector.  Use  as  many 
cards  as  needed. 

Card  1,  FORMAT  (12F6.0),  contains  the  values  of  the  elements  in  the  first  column  of  the  X 
matrix.  Use  as  many  cards  as  needed.  Repeat  this  card  for  each  additional  column  of  the  X  ma¬ 
trix. 

When  the  last  column  of  the  X  matrix  is  punched,  a  complete  problem  is  ready  to  be  processed 
by  the  computer.  The  next  card  in  the  data  file  is  Card  1,  described  above,  for  another  problem; 
therefore,  the  above  procedure  is  repeated  for  the  second  problem.  Any  number  of  problems  may  be 
processed. 

EXAMPLES  OF  DATA  INPUT 

*  *1  0  EXAMPLE  A 


3. 

22. 

39. 

54, 

65  . 

76. 

85, 

90. 

95. 

99, 

102. 

103. 

104. 

103. 

102. 

101. 

99. 

96. 

93. 

89. 

84. 

78, 

72. 

66, 

60. 

56. 

53. 

51. 

49. 

49. 

50. 

52. 

55. 

60. 

66, 

74. 

84. 

95. 

108. 

121. 

137. 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

33 
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37 

38 

39 

40 

41 

1 

4 

9 

16 

25 

36 

49 

64 

81 

100  121 

144 

169 

196 

225 

256 

289 

324 

361 

400 

441 

484  529 

576 

625 

676 

729 

784 

841 

900 

961 

1024 

1089 

1156  1225 

1296 

1369 

1*44 

1521 

1600 

1681 

1 

8 

27 

64 

125 

216 

343 

512 

729 

1000  1331 

1728 

2197 

2744 

3375 

4096 

4913 

5832 

6859 

8000 

9261 

10648  12167 

13824 

15625 

17576 

19683 

21952 

24389 

27000 

29791 

32768 

35937 

39304  42875 

46656 

50653 

44872 

5931R 

64000 

68921 

3 

28 

0  EXAMPLE 

B 

tc  LN  T  LN  T2 

269. 

180. 

150. 

120. 

90. 

75. 

66* 

60. 

54. 

48.  42. 

36. 

30, 

27. 

24. 

21. 

20* 

18. 

15.1 

14. 

12, 

10*  8. 

6. 

5, 

4. 

3. 

2.5 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1  1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1  1 

1 

1 

1 

1 

1 

0.0 

.69311 

.16321 

.60942 

.14012 

.4932? 

.72792 

.90693 

.09103 

.30323.54103 

.8177 

4.12714 

.31754 

.49984 

.76224 

.85985 

.06895 

*41615 

.57225 

.87496 

.24426.68467 

.2654 

7.62568 

.07098 

.64829 

.0119 

0. 

0.480 

1.353 

2.590 

4.580 

6.216 

7.441 

8.450 

9,55510.91112.53814.575 

17. 03318. 641 ?0. 24822. 67R23. 6 t 825. 69429. 33431. 04934, 51538. 99044. 68452. 786 
58,15065,14074,79281,214 

460  FXAMPtF  IS  TAKEN  FROM  EQUATION  31  OF  THE  TFXT. 


X2  X3  X  4 


108. 

111. 

120. 

99. 

115. 

117. 

1. 

1. 

1. 

1. 

1. 

1. 

1. 

1. 

1. 

2. 

2. 

2* 

300. 

400. 

600. 

350, 

850, 

650. 

.001 

0.001 

0.00? 

0,001 

0.00? 

0.00? 

MULTILINEAR  REGRESSION  FORTRAN  PROGRAM 

SIRFTC  SIMtIN  LIST, REF, DECK 
CMAIN 

simultaneous  linear  equations.  MAXIMUM  NUMBER  OF  ELGMENTS  IN 
COEFFICIENT  matrix  is  16,000, 

COMMON  A(16000),Yf 1000) ,0(50,50 J.YHATI 1000 ) , SS < 50 ) ,AP ( 50 ) • X ( 50 ) , 
1ST AT (50) ,I0ENT(9 J .HEAD (50) .RATIO! 50)  ,XMEAN(50) .XINVEI 50,50) ,M.N. 
2RSD0Y.FVFRE.BALM 

1  FORMAT  (2I5,F5«0.(9A6)  1 

2  FORMAT  (12A6) 

3  FORMAT  (12F6.0) 

5  format  (42H0  in  THE  GRAM-SCHMIDT  PROCESS  THFRE  WERG  .I3.30H  BINAR 
1Y  BITS  OF  ACCURACY  LOST.) 

8  FORMAT  (1H130X,(9A6)/1H08X»1HY1?X,2HY*15X.1HY12X,?HY*15X,1HY12X.2H 
1Y*15X,1HY12X,2HY*/1H  ) 

9  FORMAT  (4(1PE16.5»1°E1*«5)) 

16  FORMAT  (1H130X,(9A6)/1H051X,19HREGRESSION  ANAL YS I S/20H0  THE  MEAN 
10F  Y  IS  E1S.8/39H0  RESIDUAL  STANDARD  DEVIATION  OF  Y  IS  E15.8/29H0 
2  F  VALUE  FOR  REGRESSION  IS  F10.3/30H0  REGRESSION  COEFFICIENTS  AR 

3E/1H  1 

17  FORMAT  (8E15.6) 

18  FORMAT  (14H0  X  MEANS  ARE/1H  ) 

19  FORMAT  (30H0  INVERSE  MATRIX  ELEMENTS  ARE / 1H  ) 

?0  FORMAT  (1H  ) 

28  FORMAT  ( 1H1 30X , ( 9 A6 ) ) 

4  RFAD  1 5 »1 )M,N#SCALF» ( IDFNT ( I ) , 1*1 ,9 ) 

READ  ( 5  »2 ) (HEAD (I).I»1.M) 

READ  (5»3)(Y(|) »I*1»N) 

DO  6  J»1»M 
LV  •  (J  -  1 ) *N  +  1 
LA  »  LV  ■»  N  -  1 

6  READ  ( 3 ,3 ) ( A ( I ) . I *LV,L A ) 

SCALE  *  10,0**SCALE 
DO  13  I»1 *N 
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13  vm  -  y(I)*scale 
CALL  LINEO 

CALL  VARFAC 
LV  -  90 
DO  10  I«1*N*4 
J  »  143 

IF  (N-J)  14*19*19 

14  j  •  n 

19  IF  (50-LV)  11*11*12 

11  WRIT?  (6»8)(irFNT(K)*K>]*9) 

LV  ■  0 

12  WRITE  (6*9)((Y(K)*YmAT(K)I*K*I*J) 

10  LV  «  LV-*1 

WRIT?  (6*16MIOENT(  I)*I«1*9)*XMEAN(1 >*RSDOY»FVFRE 
LV  -  14 

DO  29  I  •  1»M*8 
J  •  1*7 

IF  CN— J»  21*22*22 

21  J  •  n 

22  WRITE  (6*17)(X(K)«K«I.JI 

29  LV  -  LVM 
WRITE  (6*18 1 

DO  30  I  ■  2 *M *8 
J  ■  147 

IF  (M-J)  23*24*24 

23  J  •  A 

24  WRITE  (6»17J(XMEAN(K)«K-I»J) 

30  LV  -  LVU 

I  •  BALM  +0,5 
WRITE  (6.5)1 
WRITE  (6.19) 

DO  29  I  ■  2.M 
IF  (LV-51 )  27*26.26 

26  LV  -  4 

WRITE  (6*2B)(ID€NT(N).N«1*9) 

WRITE  (6.19) 

27  WRITE  ( 6* 1 7 ) ( XTNVE (!*J)*J>2*t) 

WRITE  (6*20) 

29  LV  •  LV*2-*((!-2)/8) 

GO  TO  4 
END 

3IBFTC  LINEO  L I  ST .REF. DECK 
SUBROUTINE  LINEO 

COMMON  A (16000), YU 000) *0 ( 50 * 50 ) » THAT ( 1000 ) *SS ( 50 ) , AP ( 50 ) »X ( 50 ) , 
1ST*T( 50 ) ,IDENT(9)«HEAD( 50 ),RAT 10(50 ) *XMEAN( 50) »XINVE( 50*50) »M,N. 
2RSD0Y*FVERE*BALM 
DIMENSION  P(16000) 

EQUIVALENCE! A »P) 

IF  (N-l)  11*12.13 

12  X(l)  «  Y(l)/A(l) 

GO  TO  11 

13  DO  2  J-l.M 
DO  2  I-l.M 
D(I.J)  *  0,0 
IF  (I-J)  2.1*2 

1  D(I.J)  -  1.0 

2  CONTINUE 
BALM  »  0,0 
YXLO  ■  ALOG(2  *0) 

DO  24  J»2*M 
S5BFG  »  0,0 

LV  •  (J-l)*N4l 
LA  ■  LV4N-1 
DO  21  I  »  LV.LA 
21  SSBEG  »  SSBEG  ♦  P(I)*P(I) 

JM  >  J-l 
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DO  15  K-l.JK 
S  •  0*0 
R  ■  0.0 
LK  •  (K-l)*t(+i 
IK  ■  UK 

DO  *  !  »  LV.LA 
*  ■  R  +  P<I)*PILMJ 
S  •  S  ♦  P(LM)*PILM) 

♦  LK  •  LK  ♦  1 
0  •  R/S 

IF  IK-1)  20.18*20 
10  XMFANIJ)  •  0 
20  UK  •  UK 

DO  5  I  •  LV.LA 
RII)  •  Pill  -  0*PILM) 

5  LK  •  UK  ♦  1 
DO  15  I-l.K 

15  Dl  I  »  J|  ■  Dll  «  J|*»Q»D  I  T  »K ) 

SSP1K  ■  0.0 
DO  22  I  •  LV.LA 
22  SSF1N  ■  SSFIN  ♦  PI T)*P: II 
BALK  ■  BALM  ISSBEG/SSFINI 
24  COKTfKUe 

BALK-  I ALOG I BALK I /YXLO 1*0.5 
DO  7  J-l.K 
SSIJI  -  0.0 
LV  ■  IJ-ll-N+l 
LA  -  LV+K-1 
DO  A  I  -  LV.LA 
A  SSIJI  ■  SSIJ)  4  P!T)*Pin 
SSIJJ  >  SORT (SSIJI! 

DO  1A  I  -  LV.LA 
ia  phi  -  Pin/ssui 

DO  14  I-l.M 

14  DII.J1  «  Dll. J|/ SSIJ) 

API J|  *  0.0 
LK  -  0 

DO  7  I  -  LV.LA 
LK  -  LK  ♦  1 

7  APIJI  >  APIJ)  ♦  Pin*YILK) 

DO  8  I-l.M 
XII)  -  0.0 
DO  8  J-l.M 

s  xin  ■  xin-»oir.j)#Apijj 

DO  17  J-2.M 
17  Dll.J)  -  0,0 
DO  9  J-l.K 

9  STATIJ)  -  AP I J) *API J I 
DO  10  I-l.K 
YHATII)  -  0.0 
DO  10  J-l.K 
LV  •  I J— 1 ) *M4 I 

10  YHATIII  •  YHATII)  ♦  PlLV)*APlJ) 

DO  19  K-l.M 

DO  19  I-l.M 
XINVFll.r)  .00 
DO  19  J-l.M 

19  XIMVen.O  .  XINVEt  I.K)4Oll.J)*0IK.J) 

11  RETURN 
END 


AI8FIC  V ARFAC 


SUBROUTINE  VARFAC 

COMMON  *  1 16000 ) . Yt 1 000 ) .0 ( 50 . 50 ) . YhaT 1 1000 1  «SS I  AO i  ..o  i  »  vl  . 
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4  FORMAT  (1M130X.9A6) 

5  FORMAT  <?OhOSOURCF  OF  VAR t AT  TON  * 5X « I8HDE6REES  OF  FREED0M»6X»13HSUM 
IS  OF  SQUARES.16X.7HF  RATIO/1H  ) 

6  FORMAT  (7X,1A6*18X«!4.14X«E15.8*F22»3J 

7  FORMAT  IlHO,6X.8HRESIDUAl*16X«I4«14X*E13.8> 

GMEAN  »  0*0 

DO  12  1*1 *N 

12  GMEAN  ■  GMEAN*Y( 1 1 
GN  ■  N 

GMEAN  »  GMEAN /GN 
XMFANUJ  >  GMEAN 
SSTOT  ■  0.0 
DO  13  1»1 *N 

GN  «  riTt-GMFAN 

13  SSTOT  ■  SSTOT*GN*GN 
TTSTaT  «  0,0 

DO  14  J*2*M 

14  TTSTAT  >  TTSTAT+STATf J» 

SSERR  «  SSTOT-TTSTAT 
TOF  -  N-l 

ERRDE  •  TDF-fl-M 
ERRME  >  S5ERR/ERRDE 
RSOOY  «  SORT ( ERRME  j 
GN  -  M-l 

FVFRE  ■  (ERRDE*TTSTAT  T / ( GN*SSERR ) 

DO  13  I*2»M 

13  RATIOU)  ■  ST AT(  I 1 /ERRME 
WRITE  ( 6«4 1 ( I  DENT  1 1 1 » 1*1 *9 1 
WRITE  <6. 5) 

K  ■  1 

DO  16  I*2»M 

16  WRITE  l6*6JHEAO(n*K»STATIM»RATlom 
K  -  ERRDE 

WRITE  (6»T)X, SSERR 
DATA  OOOOHL/6HTOTAL  / 

HEAD! 11  ■  QOOOHl 

WRITE  (6«6>HFAO(ll*IDF. SSTOT 

RETURN 

END 
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Appendix  B 


USE  OF  THE  ORTHOGONAL  POLYNOMIAL  AND 
CURVILINEAR  REGRESSION  PROGRAM 

The  data  input  may  be  described  by  referring  to  the  curvilinear  equation 


„.£v/ 


>-  1,  2, 


where  a .  are  the  unknowns. 

Card  1,  FORMAT  (1311,  12,  313),  has  the  following  parameters: 

L  (1)  =  not  used  . 

L  (2)  =  1  to  read  in  n  ew  x  values  or  0  to  use  previous  values. 

L  (3)  =  1  to  print  out  transformation  matrix  (otherwise  0). 

L  (4)  =  1  to  print  out  orthogonal  polynomials  (otherwise  0). 

L  (5)  =  1  to  punch  out  orthogonal  polynomials  (otherwise  0). 

L  (6)  =  1  to  read  in  new  y  values  or  0  to  use  previous  values. 

L  (7)  =  1  to  print  out  regression  coefficients  for  orthogonal  polynomials 

(otherwise  0). 

L  (8)  =  1  to  print  out  conventional  regression  coefficients  (otherwise  0). 

L  (9)  =  1  to  print  out  correlation  coefficients  (otherwise  0). 

L  (10)  =  1  to  print  out  y  values  predicted  from  regression  (otherwise  0). 

L  (11)  =  not  used. 


L(8) 
L  (9) 
L  (10) 
L  (11) 
L  (12) 
L  (13) 
L  (14) 
L  (15) 
L  (16) 
L  (17) 


not  used  . 
not  used. 


L  (14)  =  m  which  is  less  than  20.  All  degrees  of  fit  from  1  to  m  are  computed. 

L  (15)  =  number  (n)  of  points.  It  is  less  than  1,000  but  must  be  greater  than  m. 

L  (16)  =  scale  factor  for  x  values  (power  of  10). 

L  (17)  =  scale  factor  for  y  values  (power  of  10). 

For  large  problems,  x  and  y  values  should  be  scaled  so  that  the  data  are  close  to  unity. 
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Card  2,  FORMAT  72H  .  .  .  ),  identifies  the  problem  with  any  legal  Hollerith  characters.  The 
character  in  the  first  column  enters  the  carriage  control  channel. 

Card  3,  FORMAT  (12F6.0),  contains  the  *  values.  Use  as  many  cards  as  needed.  This  card 
must  not  be  in  the  data  file  if  L  (2)  =  0,  in  which  case  the  s  values  from  the  previous  problem 
are  used.  (L  (2)  =  1  for  the  first  problem.) 


Card  4,  FORMAT  (12F6.0),  contains  the  y  values.  Use  as  many  cards  as  needed.  This  card 
must  not  be  in  the  data  file  if  L  (6)  =  0,  ii*  which  case  the  y  values  from  the  previous  problem 
are  used.  (L  (6)  =  1  for  the  first  problem.) 


With  the  completion  of  Card  4,  the  next  card  in  the  data  file  is  Card  1,  described  above,  for 
the  next  problem;  thus/ the  above  process  is  repeated  for  another  problem.  Any  number  of  prob¬ 
lems  may  be  processed. 


The  program  is  coded  in  single  precision,  and  thus  large  differencing  errors  may  appear  for 
high  degrees  of  fit. 


EXAMPLES  OF  DATA  INPUT 


1  1  111  114-6-4 

1  EXAMPLE  ONE 

3002  2915  2832  2 754 

27993  25224  22989  20569 

1  111  1  4-6-4 

1  EXAMPLE  TWO 
28195  25623  23201  20792 
1  111  14-6-4 

1  EXAMPLE  THREE 

3007  2915  2826  2743 

111  1  45 

1  EXAMPLE  FOUR 

-2-1  0  1  2 

1  2  3  4  5 

1  1  111  4  10  -5  -3 

1  example  five 

2525  2828  3147  3586  4198  4456  5757  6688  7662  8704 

1000  1099  1144  1188  1228  1253  1304  1339  1374  1408 
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ORTHOGONAL  POLYNOMIAL  AND  CURVILINEAR  REGRESSION  FORTRAN  PROGRAM 


SINFTC  CURVLN  L I  ST  .REF.OECIC 

COWON  P(999»20) »M,N,L(17).Y(999)*  B(20).  SS(20|.  X8(20,20).  TOL 

1  FORMAT(13ll»I2*3I3»F3»6) 

2  FORMAT I 54H1BAD  CONTROL  CARO  OP  VISCOUNT  OF  INPUT  OUlT  THIS  JOB.) 
10  READ  (3«l)(L(l),t  ■  1*17) 

9  M  •  L(14) 

N  •  L(13>  , 

IF  (N-l-M)  1111.13.13 

13  DO  14  K  •  1.13 

IF  (LOCI  -1)  14.14.1111 

14  CONTINUE 

IF  (M-19)  13.13.1111 
13  CALL  ID 

IF  (L(2) )  17.17.16 

16  CALL  XFILE 

17  IF  (L ( 3  ) )  19.19.18 

18  CALL  OETMAT 

19  IF  (L(4))21.21.20 

20  CALL  WRITEP 

21  IF  (L(5))26.26.24 
24  CALL  PUNCHP 

26  IF  (L (6 )  )  28.28.27 

27  CALL  VFILE 

28  IF(L(7))30. 30.29 
2°  CALL  OUT 8 

30  IF (L ( 8 ) ) 32 .32 .31 

31  CALL  OUTA 

32  IF(L(9))34.3*.33 

33  CALL  OUTR 
CALL  OUTFR 

34  IF(L(10)  110. 10,33 

35  CALL  OUT 
GO  TO  10 

1111  WRITE  16.2) 

CALL  EXIT 
END 

SIBFTC  OUTFR  LIST, REF. DFCX 
SUBROUTINE  OUTFR 

COMMON  PI999.20) ,M.N»L( 17) »Y(999) .  B  <  2  0 ) •  SS(20).  XBI  20.20),  TOL 
COMMON  SSY 

DIMENSION  FA ( 20  )  ♦  FR(20),  0 ( 20 ) »R ( 20) . S ( 20 ) 

1  FORMAT  ( 51H0F— RATIOS  (INDIVIDUAL)  D.F.(NUV)»1»  D.F. (DEN ) »N-1- J) 

2  FORMAT (65H0F-RATIOS  (ACCUMULATIVE)  D.F.(NUM)«J,  D.F. ( DEN ) =N-1-J .  J 

1»0,1 «•••»,  12 ) 

3  FORMAT  (1H  ) 

4  FORMAT! 1H0.10F11 .3) 

7  FORMAT  (25H0RESIDUAL  VARIATION  IN  Y.) 

8  FORMAT  (1H0.7E15.6) 

9  FORMAT  (36HORESIDUAL  STANDARD  DEVIATION  OF  FIT,! 

00  *  N-2 

MP  »  M+l 

DO  5  J  «  2.MP 
3  0(J)  *  B(J)**2 
D  «  SSY 
A  *  0.0 

DO  6  J  *  2.MP 
S  «  N  -  J 

T  «  J  -  1 
00  «  S/T 
A  *  A  ♦  Q(J) 

D  ■  D  -  0(JI 
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R(J)  •  o 

faij)  ■  o(J)/o  *oo 
6  FB(J)  •  A/D*QQ 
R(l)  •  SSY 
00  11  J  «  l.MP 
X  -  N-J 

11  S(J)  •  SQRT(R(J)/X) 

WRITE  (6.3) 

WRITE  ( 6*1 ) 

WRITE  (6.4)(FA(J).  J  ■  l.MP) 

WRITE  (6.3) 

WRITE  (6.2 )M 

WRITE  <6.4)(FB(J>.  J  »  l.MP) 

WRITE  (6.3) 

WRITE  (6*7) 

WRITE  (6*6 ) (R ( J) .  J  -  l.MP) 

WRITE  (6.3) 

WRITE  (6.9) 

WRITE  (6.8)(S(J).  J  ■  l.MP) 

RETURN 

ENO 

SIBFTC  OUTR  L 1ST .REE. DECK 

SUBROUTINE  OUTR 

COMMON  P(9R9*20) »M.N.L( 17) »Y(999) .  8(20).  SS«20).  XB(20.20)»  TOL 
COMMON  SSY 

DIMENSION  R(20).V(20) 

1  FORMAT  ( 26HOCORRELAT  TON  COEFFICIENTS.) 

2  FORMAT  (5°hOY  VARIATION  ACCOUNTED  FOR  3Y  J  TH  POLYNOMIAL  SET  ) 

3  FORMAT  (1H  ) 

*  FORMATI1H  .F9.3.11F10.5  ) 

3  FORMATI1H  .7E13.6I 
WRITE  (6.1) 

WRITE  (6.3) 

0  ■  1.0/SORT (SSY) 

MP  •  M  ♦  1 
DO  6  J  »  2.MP 
V(J)  »  B ( J) **2 

6  RiJ)  ■  B ( J ) *0 

WRITE  (6.4 ) (R I J) . J  »  1*  MP) 

WRITE  (6.3) 

WRITE  (6.2) 

WRITE  (6.3) 

WRITE  (6.3)(V(J).J  ■  l.MP) 

RETURN 

ENO 

SIRFTC  ID  LIST. REF. DECK 

SUBROUTINE  ID 

COMMON  P ( 999 . 20) .M.N *L 1 17 ) *Y ( 999 ) ♦  8(20).  SS(20)»  XB(20.20).  TOL 
COMMON  SSY 

1  FORMAT (63H10RTH060NAL  POLYNOMIAL  ANO  CURVILINEAR  REGRESSION  ROUTl 
1NE  131-22. 112. 1111) 

2  FORMAT  ( 72H 

1  ) 

5  FORMAT  (1H  ) 

3  FORMAT ( 13HO0FGRFF  FIT  -,[3»5H»  N  -.14) 

4  RE AO  (3.2) 

WRITE  (6. 1ML(I). 1-1.12) 

WRITE  (6.2) 

WRITE  ( 6 .3 ) M.N 
WRITE  (6.5) 

RETURN 

ENO 

SIRFTC  DFTmaT  LIST. REF. OFCX 
SUBROUTINE  DETMAT 

COMMON  P(999«20) .M.N »L ( 1 7 ) . Y ( 999 ) ♦  B(20)»  SS<20).  XB(20.20).  TOL 
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COI#K)N  SSY 

1  FORMAT  (?4H0MATRIX  OF  DETERMINATION) 

2  FORMAT  (1H  ) 

3  FORMAT  (1H  .7E1S.6) 

6  FORMAT  (  1  HI  1 

A  WRITE  (6*1) 

WRITE  (6*2 ) 

MP  •  M  ♦  1 

DO  3  J  »  l.MP 

WRITE  (6,3>(X8(K,J)  .K  ■  1*J) 

3  WRITE  (6*2) 

WRITE  (6*6) 

RETURN 

END 

SIRFTC  WRITEP  LIST, REF. DECK 
SUBROUTINE  WRITEP 

COMMON  P(999»20),M.N,L<17)»Y(999).  8(20),  SSI20)*  X8(20*20)»  TOL 
COMMON  SSY 

1  FORMAT  ( 23H0ORTHOGONAL  POLYNOMIALS) 

3  FORMAT  (5HOJ  >  *1?) 

A  FORMAT  (1H  * E20# 6 *6E1 5 *6 ) 

3  WRITE  (6*1) 

MP  ■  M  ♦  1 

DO  6  J  »  2»MP 
K  -  J  -  1 

2  WRITE  (6*3 )K 

6  WRITE  (6*A)(P(I.J)  *  I  -  1 »N ) 

RETURN 

END 

SIBFTC  OUTB  L I  ST , REF. DECK 

SUBROUTINE  OUTB 

COMMON  P(999*?0) ,m,n.L( 17)*Y(999) ,  8<20).  SS(20)»  XB(20»20),  TOL 
COMMON  SSY 

1  FORMAT  ( 52HOREGRESS ION  COEFFICIENTS  FOR  ORTHOGONAL  POLYNOMIALS  ) 

2  FORMAT  (1H  ) 

4  FORMAT  (1H  *7E13«6) 

3  WRITE  (6*1) 

WRITE  (6*2) 

MP  ■  M  ♦  1 

WRITE  (6.A) (6( J) *  J  ■  l.MP) 

WRITE  (6.2) 

RETURN 

END 

JIBFTC  XFILE  LIST. REF, DECK 
SUBROUTINE  xfile 

COMMON  p(999.20) ,m,n,L(17).Y(999).  S(20).  SS( 20),  XB(20.20).  TOL 
COMMON  SSY 
DIMENSION  R2I20) 

1  FORMAT  (12F6.0) 

2  READ  (3.1 ) (P ( I »2) .  1*1. N) 

MP  *  M  +  1 


3 

32 

33 
36 


5 

A 


DO  3  I  ■  1 ;N 
P(I.l)  *  1.0 
IF  (3-MP)  32.32.36 

DO  33  J  »  3 .MP 
DO  33  I  *  1 »N 
P ( I • J )  «  P(l*2>**( J-l) 

DO  A  J  •  l.MP 

DO  A  I  ■  l.MP 

X8II.J)  «  0.0 

IF  (I-J)  4,5.4 

XB(I.J)  »  10,0m(L(161*(1-J)) 

CONTINUE 

DO  7  J  *  2.MP 

JM  •  J-l 
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00  7  K  m  1,jm 
S  m  Ot 0 
R  m  0.0 
00  6  I  •  l.N 
«  •  R  4  P(( 

A  s  •  S  4  P(  r,K J#*2 

c  ■  n/s 

DO  8  I  •  l.MP 
•  :  »•'->>  - 
7  SV;  plJ,J>  m  c»p/r»tn 

00  9  J  •  1 ,MP 
SSI  J)  m  0.0 

00  9  ?  •  ltK 

*  nni!  m,SSlJt  * 

00  n  j  m  l.MP 

5  •  UO/SORTfSSiJtt 

OO  10  !  •  1.MP 

10  XSd.J)  m  XBll  ,J)#S 

00  H  I  .  1,N 

11  b«i.j)  •  Ptujfs 
return 

END 

«-!ST  .REE, DECK 
SWWlffHE  YFfLE 

COMwJI  ss?9<’,?0,*M*N*‘-n7».Y(999>.  B<20>,  SSI20)  ,  X8(20.20).  Tot 
DIMENSION  YY ( 999 )  01 

j  "«*■«» .  sUM  of  ,«,««  y  . 

3  FORMAT  (12E6.0) 

REXO  rs.3)  (Y(n  »I«1»N) 
o  •  io,o**Ln?) 

00  4  r  «  l.N 
Y(l»  mo*  Yt II 

4  mn  >  Ym 

HP  •  M  4  1 
00  9  J  •  l.MP 
BIJ)  .  o 
00  9  I  m  l.N 

5  Sojj 

’  Ur".  o  rT'"  ‘  *<■»«■«•■» 

YB  «  B«l)*XBIl,i  ) 

00  o  r  •  i.n 

4  WRITE  fJl7)SSY,n  *  Y8’  *#2 
WRITE  (8.2) 
return 
END 

,,"'TeJ5SOT„i'5;r*"c' 

iss:  sY’-’0’ . iT..r.w„  „„„  I0C 

DIMENSION  4(20.20  ) 

i  to™”  iiiToM^nV.  :nrs'°*  c°^na^^*i 

?  ass  jar™-* 

write  (6,n 

WP  *  M4l 

4JI1*1)  «  Bll)  *XB  (1,1) 

00  6  J  m  2.WP 

00  4  X  »  l.j 
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4  A(J»X)  •  A(J-l.X)  ♦  XB(X*J)*eU> 

6  CONTINUE 

00  5  J  »  1»MP 
XM  ■  J-l 
WRITE  (6«2)XM 

3  WRITE  (6*3  HA  I  J*K)  ,  K  ■  1*J) 

WRITE  (6*7) 

RETURN 

ENO 

SIBFTC  OUT  LIST, REF. DECK 
SUBROUTINE  OUT 
DIMENSION  T( 10 ) 

COMMON  P(99R*20) *M.N.L<17)»Y(999) t  B(20>»  SS(20)*  XB(20*20)*  TOL 
COMMON  SSY 

1  FORMAT  (39H1Y  AND  PREDICTED  Y«S  FROM  M  DEGREE  FIT.) 

2  FORMAT  (1H  1 

3  FORMAT  (1H  ♦F19«0«F10,O.F20,0,F10I0,F20,0*F10,0»F20.0»F10.0) 

A  FORMAT  (13HO{NOT  SCALED)) 

WRITE  (6.1) 

WRITE  (6*4 ) 

WRITE  (6*2) 

0  •  1.0/10.0**L(17) 

MP  •  M-.1 
XL  *>  6 

00  6  I  ■  1*N*4 
IF  (N-I)  9.8*9 

8  XL  *  2 

9  T  ( 1 )  »  Y  ( I )  *0 
T( 2 )  ■  0.0 

IF  IN-I-1)  11. 10,11 

10  XL  •  4 

11  T ( 3 )  *  Y(  I  +1 )  *0 
T{ 4 )  «  0.0 

IF  (N— 1-2)  13.12,13 

12  XL  -  6 

13  T(  5 )  *  Y(I-*2)*0 
T ( 6 1  «  0.0 
T(7)  •  Y( 1*3) »0 
T ( 8 )  ■  0.0 

DO  5  J  «  1 .MP 
T ( 2 )  «  T (2 )  +  P< I.J)*8< J)*0 

T ( 4 )  «  T(4 )  ♦  P( I+1*J)*8(J)*Q 
T(6)  »  T ( 6 )  4  P( I42«J)*B( J)*0 
3  T  ( 8  )  «  T  (8  )  P(  1+3,J)*6(  J)*0 
6  WRIT?  (6*3)(T(K).  X  »  1,8) 

RETURN 

ENO 

SIBFTC  PUNCHP  LIST .REF.DECX 
SUBROUTINE  PUNCHP 

COMMON  P(999,20) ,M,N.L( 17) ,Y(999) ,  B(20)»  SSI20),  XB(20.20),  TOL 
COMMON  SSY 

1  FORMAT  (23HOORTHOGONAL  POLYNOMIALS) 

3  FORMAT  (3H0J  *  *12) 

4  FORMAT  (1H  ,E 1«.6»4E19,6 ) 

5  WRITE  (3*1) 

MP  ■  M  ♦  1 

DO  6  J  «  1»MP 
K  •  J  -  1 

2  PUNCH  3,  K 

6  PUNCH  4,  IP(!*J)»  I  *  1 »N ) 

RETURN 

END 
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Appendix  C 

USE  OF  THE  FACTORIAL  ANALYSIS  PROGRAM 

The  data  input  may  be  divided  into  two  parts.  The  first  part  is  the  orthonormal  tensors,  and 
the  second  part  is  the  measured  data  points  which  are  to  be  arranged  in  a  factorial  design.  Using 
the  notation  of  the  text,  the  orthonormal  tensors  are  designated  as 
2  <»  <9 

fl‘  1  <s  <i 

'  1  <  /  <  i 

The  index  t  designates  a  tensor  which  is  congruent  to  a  square  matrix  of  order  t.  The  index  a 
denotes  a  normal  vector  of  the  tth  tensor,  and  the  index  /  denotes  a  component  of  the  sth  normal 
vector.  For  an  tth  tensor,  there  are  i  normal  vectors  and  i  components  for  each  normal  vector. 
Thus  there  are  t2  elements. 

Card  1,  FORMAT  (1 1 ),  gives  the  order,  t,  of  the  tensor  which  appears  on  the  next  card. 

Card  2,  FORMAT  (5X,  5E  15.7),  contains  the  »  values  of  the  first  normal  vector  of  the  tth 
tensor.  If  i  is  greater  than  5,  use  another  card  with  the  same  format.  Repeat  this  card  for  each 
additional  normal  vector  of  the  tth  tensor. 

When  the  ith  normal  vector  is  punched,  repeat  the  above  cards  for  another  tensor.  The  process 
is  terminated  when  t  is  equal  to  zero.  These  tensors  (orthogonal  polynomials)  may  be  obtained 
in  the  proper  format  from  the  program  in  Appendix  B. 

The  second  part  of  the  data  input  is  for  the  measured  data  points  of  a  factorial  design. 

Card  1,  FORMAT  (215,  F5.0,  57H  .  .  .)  has  the  following  parameters: 

NO)  =  number  of  factors  in  the  design.  It  must  be  less  than  15.  Replication  is  not 
considered  to  be  a  factor. 

N  (2)  =  number  of  data  points  in  the  design.  It  must  be  less  than  21,001. 

SF  =  scale  factor  for  the  data  points  (power  of  10). 

ID  =  identification  of  problem  with  any  legal  Hollerith  characters. 

Card  2,  FORMAT  (4X,  1514),  has  the  following  parameters: 

I.  (1)  =  number  of  replications. 

I,  (2)  =  number  of  levels  for  first  factor. 

L  (3)  =  number  of  levels  for  second  factor. 

L  (15)  =  number  of  levels  for  fourteenth  factor. 

If  a  minus  sign  is  punched  in  front  of  the  number  of  levels,  that  factor  is  considered  to  be 
quantitative.  Otherwise,  the  factor  is  considered  to  be  qualitative.  Each  distinct  number  punched, 
for  the  number  of  levels  of  a  factor,  presupposes  that  an  orthogonal  tensor  of  that  order  has  been 
included  in  the  first  part  of  the  data  file. 
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Card  3,  FORMAT  (12F6.0),  contains  the  measured  data  points  of  the  factorial  design.  Use 
as  many  cards  as  needed.  The  data  points  must  be  ordered  by  the  level  indices  of  the  replicate 
and  of  the  factors,  according  to  the  hybrid-base  number  system  (see  text),  where  the  highest- 
order  index  is  for  replication,  the  next  order  index  is  for  the  first  factor,  etc.,  until  the  lowest 
order  index  is  for  the  last  factor  of  the  design. 

When  all  the  data  points  have  been  punched,  one  factorial  design  has  been  completed.  Any 
number  of  factorial  designs  may  be  processed  by  repeating  the  above  cards  for  the  second  part  of 
the  data  input.  The  first  part  of  the  data  input  is  not  repeated  for  each  factorial  design.  If  dif¬ 
ferent  tensors  are  to  be  used  for  a  design,  a  new  data  run  must  be  made  on  the  computer  with  both 
parts  of  the  data  input  present.  A  limitation  on  the  use  of  the  program  is  that  two  factors  having 
the  same  number  of  levels  cannot  be  analyzed  with  two  different  tensors. 


EXAMPLES  OF  DATA  INPUT  !  , 

2 

7*07106783-01  7.07106783-01 

—7 *071 0678 P-01  7*07106783-01 

3 

5*77350273-01  5.77350273-01 

-7.07106783-01  0, 

4.0824829F-01  -6,16496563-01 

4 

5*00000003-01  3*0000000^-01 

—6*70870393-01  -2*73606803-01 
3  *  OOOOOOOf-0 1  -5.00090003-01 
-2*23606603-01  6.7082039F-01 

3 

4*4721 360p-0 1  4.4721360E-01 

-6.3245353F-01  -3, 16727773-0} 
3.3432248F-01  -7 ,6726 1 24 3-01 
— 3 • 1 622777F-0 1  6* 32435 63F-0 1 

1 « 1952286F-0 1  -4,78091443-01 

6 

4*08748 79 F-01  4,06?4879F-01 

4.0324879F-01 

— 3*976 1470 F— 01  -3.5856e583-0l 
5  *97614  30F-0 1 

5 *45544 73 3—0 1  -1  *091 08953-0 1 
5*45544733-0 1 

-3.7767800F-01  5. 2 1 749193-01 

3*72678003-01 

1 .8398724F-01  -3,6694671 F-01 
1.8898224F-01 

-6* 7994079 F-0 ?  3«1497039F-01 

6 *2994 079 F-0 7 
7 

3*77964473-01  3.77964473-01 

3*7796447E-01  3 . 7796447E-01 

—3  *6694671 F-0 1  -3  »  77964 4 7F-0 1 
3«7796447F-0l  5,66946713-01 

5*4554473F— 0 1  0, 

0,  5«4534473F-01 

-4.0824829F-01  4,08248293-01 

-4«0824829F-0 1  4.0824829F-01 

2 *41 74689F-0 ]  -5 *6 40 76 07 F-01 
-3 , 640 76 07 F-0 1  2,41 74689F-01 

-1 *091 0895 F-01  4,3643578F-0] 

-4, 36* 3 378 F-0 1  1 .0910895F-01 

3,?897383F-02  -1  * 97783 3 1 F-01 
— 1 *9738351 F-0 1  3,2897585F-02 


5«7735027F-01 
7*071 0678F-01 
4.0824829E-01 

5*09000003-01  3.0000000E-01 

7*23606803-01  6 . 7082039E-01 

-5.0009000F-01  5  *  00000003-01 

-6.7082039E-01  2.2360680E-01 

4,47213603-01  4.4721360E-01  4. 4721 360E-0 1 

0,  3.1622777E-01  6. 3245553E-01 

-5.3452248E-01  -7 . 67261 24E-01  5. 34522483-01 

O*  -6.3243353E-01  3. 1622777E-01 

7,171 3717F-01  -4.7809144E-01  1 . 1 952286E-0 1 

4.0824829F-01  4.0824829E-01  4.0824829E-01 

-1*19322863-01  1. 1932286P-01  3,58568583-01 

-4.3643578E-01  -4 . 3643578E-0 1  -1 ,0910895E-01 

2  «  98 1 4240E-01  -2.9814240E-01  -5. 2 17491 9E-01 

3*7796447F-01  3.7796447F-01  -5.6694671E-01 

-6*29940793-01  6.2994079E-01  -3. 1497039E-0 1 


5.7796447E-01  3.7796447E-01  3. 7796447E-0 1 

-1.3898224E-01  0.  1.8898224E-01 

■3*27?2684E-01  -4. 3643578E-01  -3.2732684E-01 
4.0824829E-01  0.  -4.0824829E-01 

8*058?296E-02  4,83493783-01  8.0382296E-02 

•5*45544733-01  0,  5.45544733-01 

4*93463773-01  -6*57971693-01  4, 9346377E-0 1 
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9 


0 


3.5355339E-01 
3.3355339E-01 
-5 .4006 172 E-Ol 
2* 3145502E-01 
5.4006172E-01 
-2.3145502E-01 
-4.3082022E-01 
-4.3082022F-01 
•2.8203804F-01 
-1.2087344F-01 
-1.4978617E-01 
3.6376642E-01 
6.154S745E-02 

5.5391171E-01 

-1.7069719E-02 

3.5846409E-01 

3.3333333E-01 
•3«3333333E-01 
“5  «1 639778 E-Ol 
1.2909944E-01 
5.3181602E-01 
-3.2288830E-01 
-4,44*4421 E-Ol 
-2.8603878F-01 
3. 1289311 E-Ol 
2.01 14937E— 0 1 
-1.8490007F-01 
4.1602919E-01 
8.9893319E-02 
2 .2473329F-0  2 
-3.4139437E-02 
-4.7795212E-01 
8 .81 476*8 F-0 3 
-4.9362683E-01 


3.9399339E-01 
3.5399339E-01 
-3.8579837E-01 
3.8979837E-01 
7.7191675E-02 
7.7191679E-02 
3.0772873E-01 
-3 .0772 873 E-Ol 
-9. 23 79493 F-0 1 
-9.2379493F-01 
4.921 9497E-01 
-4.9219457E-01 
-3. 0772 8 73 E-Ol 

-3.0772873E-01 

1.1948803E-01 

-1.1948803E-01 

3.3333333E-01 
3.3333333E-01 
-3. 8 7298 33 E-Ol 
2.9819889E-01 
1 .3299*01 E-Ol 
-1.9194744E-01 
2.2247460E-01 
— 4.1316712E-01 
-4.6933966E-01 
-2.4584439E-01 
5.0647518E-01 
1 . 8490007 E— 01 
-3. 820*6 59 E-Ol 
4.94*1 323E-01 
2.0483662F-01 
4.7795212F-01 
-7.0518 118E-0 2 
2.4681 341E-01 


3.5359339E-01 

3.5355339E-01 

-2.3145502E-01 

5.4006172E-01 

-2.3145502E-01 

5.4006172E-01 

4.3082022F-01 

4.3082022F-01 

-1.2087344E-01 

2.8203804E-01 

-3.6376642E-01 

1.4978617E-01 

3.3391171E-01 

6.1545745E-02 
-3. 5846*09 E-Ol 
1.7069719E-02 

3.3333333F-01 
3.3333333E-01 
-2. 58 1 98 8 9 E— 01 
3.8729833E-01 
-1.5 194744 E-Ol 
1.3295401E-01 
4.1316712E-01 
—2. 224746 0F-01 
— 2.4564459E— 01 
-4 .69 3 3966 E“01 
-1.8490007F-01 
-5«0847518E-01 
4.9441323E-01 
-3.8204659E-01 
-4.7799212F-01 
- 2. 048 3662 E -01 
2.4681341E-01 
-7.0518118E-02 


3.5355339E-01 
-7.7151675E-02 
-3.8575837E-01 
1 . 8463724E-01 
3.6262033E-01 
-3.2097037E-01 
-3.0772873E-01 


5.9744015E-01 


3. 3333333F-01 
3.3333333E-01 
-1 »  2 90994* E— 01 
5.1639778E-01 
-3.2288830E-01 
5  »  3181602E-0 1 
2.8603878E-01 
4 .4494921 E-Ol 
2.0114557E-01 
3. 128931  IE-01 
-4. ) 60251 5E-01 
1 . 6490007E-0 1 
2.2473329E-02 
8.9893315E-02 
4. 77952 1 2F-01 
3.4139437E-02 
-4.9362683E-01 
8.8147648E-03 


3.5393339E-01 
7.7131675E-02 
-3.8975837E-01 
-1.6463724E-01 
3.6262033E-01 
3.2097037E-01 
-3 .0772873 E-Ol 

-5. 974401 5E-01 

3.3333333E-01 

0. 

-3. 7986859E-0 1 

0. 

4.0229114E-01 

0. 

-4.4946657E— 01 

0. 

6.1703353E-01 


4 
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-3  EXAMPLE 

NUMBER  ONE 

1  2 

2  ■ 

-3  -4 

2793 

2007 

1959 

3010 

2913 

2237 

2649 

2225 

1990 

2582 

23*2 

2210 

2438 

2312 

2086 

2*05 

2766 

2419 

2339 

2225 

1944 

2838 

2579 

2003 

4 

189 

0  EXAMPLE 

number  two 

1  -3 

7-3  3 

1.209 

1.206 

1.201 

1.203 

1.205 

1.204 

1.200 

1.203 

1.201 

1.200 

1.206 

1.203 

1.209 

1.210 

1.212 

1.150 

1.148 

1.150 

1.149 

1.138 

1.139 

1.140 

1.139 

1.139 

1.163 

1.165 

1*161 

1.165 

1*162 

1.163 

1.049 

1.052 

1.053 

1.050 

1.043 

1.043 

1.030 

1.0*9 

1.050 

1.050 

1.050 

1.0*9 

1.315 

1.319 

1.319 

1.305 

1.339 

1.320 

1.330 

1.328 

1.325 

1,320 

1.319 

1.322 

1.390 

1.383 

1.383 

1,385 

1,384 

1*388 

1.399 

1.395 

1.397 

1.395 

1.398 

1.395 

1.395 

1.392 

1.391 

1,180 

1.185 

1.180 

1*169 

1.139 

1.161 

1,165 

1.163 

1,161 

1.173 

1.178 

1.178 

1.180 

1.179 

1.179 

2.000 

1.994 

1.999 

2.020 

2.012 

2.015 

2.020 

2.021 

2.023 

2.020 

2.019 

2.025 

1973 

2563 

2666 

2160 

1803 

2520 

2105 

2764 

2503 

2215 

2064 

2843 

2325 

2644 

2608 

2203 

2048 

2544 

2218 

2*94 

2655 

2230 

2084 

2712 

1.210 

1.208 

1.210 

1.200 

1.199 

1.200 

1.205 

1.206 

1.205 

1.205 

1.205 

1.210 

1.150 

1.151 

1.155 

1.155 

1.153 

1.156 

1.140 

1.140 

1,140 

1.160 

1.159 

1.159 

1.050 

1,056 

1,050 

1.055 

1.053 

1.059 

1,050 

1.042 

1,041 

1.035 

1.043 

1.0*5 

1.035 

1.050 

1.050 

1.320 

1.318 

1.319 

1.330 

1.326 

1,314 

1.320 

1.322 

1.325 

1.315 

1.312 

1.320 

1.320 

1.312 

1.319 

1.390 

1,385 

1.390 

1.395 

1.397 

1.393 

1.400 

1.395 

1,391 

1.395 

1.395 

1.391 

1.180 

1.160 

1,180 

1*185 

1.182 

1.182 

1.165 

1.161 

1.163 

1.180 

1.175 

1.175 

1.990 

1,988 

1.990 

1.990 

1.990 

1.993 

2.015 

2.016 

2.016 

2,020 

2.015 

2.016 

2,025 

2,022 

2.025 
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FACTORIAL  ANALYSIS  FORTRAN  PROGRAM 


•T8FTC  MAIN 
CMAIN 


LIST.REF.DECK 


ANALYSIS  OF  VANfANCF 

1  FORMAT  (2I5.E5.0 ,9A6« A3 ) 

2  FORMAT  I AHO  ,15141 
1  FORMAT  (12F6.01 

4  FORMAT  IflOnO  Brp  •  *  yp  *«  v. 

x13  X1A  SUM  OF  SOUARES/1H  T  ”  6  7  X*  X9  X1°  X11  X 

5  FORMAT  UH  .15I4.E16.81 
15  FORMAT  1 1M1 »30X » 10A6 ) 


'  '  ii:  _  15 *o I 

15  FO#WAT  ( 1H1 #30X* 10A6 ) 

Ki?w  atpass.'s  :;!:}*•  •>»•«»  »*«  -»■«>  of  m, 


call  rolt  ' 

IEi0N^!5*1,M*N*SCALP'<Me*0<n,,“1*9,*Hf*B<l0> 

RFAO  (5»2)(LEV(F) ,  K  .  i,*> 


RFAO  I5.2HLEV(K),  k  .  l,N) 
00  23  K«1,M 
TF  fLEV(*n21  .13.22 

21  LINIK)  a  l 

LFVIkj  ,  -LEVI*) 

60  TO  23 

13  LEV  (F)  .  1 

22  LIN  (K )  a  0 

23  CONTINUE 

14  IF  (LEVUI-l)  18.18.20 

18  MF  ■  2 
GO  TO  19 

20  MF  •  1 

19  SC  »  SCALE 
LUMP ( M )  *  1 
00  7  I  a  2.M 
J  -  M-I-M 

7  LUMP! J)  a  LUMP( J>1 )*LFVf J+l  ) 
READ  I5.3IIDIJI.  I  a  i,N) 
SCALE  a  10,0#*$CALE 
DO  9  I  »  1,N 
9  Dill  a  D(I)*SCALE 
00  17  FM.15 
17  FOG(K)  »  0 

WRITE  (6.15HHEA0I  I).I  *  1,10) 

K  a  M-l 

WRITF  (6.161K.N.5C 
CALL  MEANS 
CALL  ANLSIS 
00  12  I  a 1  ,N 

12  01 1 1  a  01 1 )**2 

LEVEL  «  1 
LV  a  50 
00  8  I  a  l.N 
IF  I 50-LVI  10,10,24 
10  WRITfc  ( 6  *  1 ! > (MEAOIJ) » J  -  1,101 
WRITE  (6.161F.N.SC 
WRITE  I6.2)(LEV(K),  K  •  1»M> 
WRITE  (6,4) 

LV  ao 

24  WRITE  <6.5l(KOG<K)»  K  a  1H5), 
CALL  COUNT  (16) 

8  LV  a  LV+1 
CALL  ANOVA 
GO  TO  6 
END 


»  1*15).  DU) 
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trsrrc  POLY  LIST .REF.DECK 

subroutine  poly 

. . 


1  FORMAT  fH) 

2  FORMAT  (1H  • 

3  READ  <5. 1)N 
IE  (N)  6*6.4 

A  DO  5  J  •  1*N 

5  READ  f  5  *?) (P (N* J »  U  > 
60  TO  3 

DATA  OOOOHL/6HREP  / 

6  HOLER !  1 )  ■  OOOOHL 

DATA  0001HL/6HX1  / 

HOLER !  2 )  ■  0001HL 
OATA  Q002HL/6HX2  / 

H0LER!3)  »  0002HL 
DATA  0003HL/6HX3  / 

HOLER!*)  b  Q003HL 
DATA  Q00AHL/6HX*  / 

HOLER !  5 )  a  0004HL 
OATA  0005ML/6HX3  / 

HOLER !6)  ■  Q005HL 
DATA  0006ML/6HX6  / 

HOLER (7)  a  Q006HL 
DATA  0007HL/6HX7  / 

H0LER<8)  ■  0007HL 
DATA  G008HL/6HX8  / 

HOLER ( 9 )  a  Q006HL 
OAT*  OOO9HL/6HX9  / 

HOLER! 10)  a  Q009HL 
DATA  0010HL/6HX10  / 

HOLER! 11)  »  OOlOHL 
DATA  0011HL/6HXH  / 

H0LERI12)  a  00 UHL 
DATA  0012HL/6HX12  / 

HOLER  113)  a  QOuhl 
DATA  001 3HL/6HX1 3  / 

H0LER11A)  a  0013HL 
DATA  001AHL/6HX1A  / 

HOLER  1 13)  a  001 AHL 
DATA  0015HL/6HLIN  / 

H0LEPU6)  »  Q015HL 
DATA  Q0l6ML/«TtoUA  / 

HOLER! 17)  a  0016HL 
DATA  0017HL/6HCUB  / 

HOLER! 18  I  a  Q017HL 
DATA  0018HL/6H.  / 

HOLER ! 1 9 )  a  Q018HL 
DATA  0019HL/6HI  / 

HOLERI20)  a  0019HL 

DATA  0020HL/6H).  / 

HOLER  I  21)  s  0020HL 
DATA  0021HL/6H  ,/ 

HOLER! 22)  a  002] HL 
RETURN 
END 

SIBFTC  ANLSIS  LIST  .REF.DECK 
SUBROUTINE  ANLSIS 


E19.6.  4E13.6I 


I  •  l.N) 


1 


00  1  K  B  1,13 
K06 ! K )  a  0,0 
LEVEL  ■  LEV! J ) 

00  3  I  a  l.N, LEVEL 
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CAUL  COUNT  (J) 

00  2  K  *  1 .LEVEL 
IN  •  INDEX OO 

2  VECTOR(K)  «  D(IN) 

DO  3  K  •  1. LEVEL 
IN  •  INDEX(K) 

DC  IN »  •  0.0 

DO  3  L  »  1. LEVEL 

3  OfINI  «  0(  IN)  ♦  VECTOR! L)*P( LEVEL »K,L> 

RETURN 

END 

SIBFTc  COUNT  LIST. REP. DECK 

SUBROUTINE  COUNT  (J) 

COMMON  P (9.9.9) .INDEX ( 15 ) .KOGI 15) .LEVI  15) .LUMP! 15 ) .01 21000 J .LIN ( 15 
1 ) .VECTOR (9). HEAD (10), HOLER (22 ) .LEVEL. M.N.MF 

1  INDEX! 1)  ■  1 
DO  2  I  •  1 .M 

2  INDEX ( 1 )  •  INDEX ( 1 )  +  LUMP ( 1 ) *KOG (ft 
DO  3  I  *  2 .LEVEL 

3  INDEX! I)  ■  INDEX ( I -1)  ♦  LUMP(J) 

DO  5  I  *  l.M 

L  «  M-IM 
IP  (L-J)  A. 3. A 
A  KOG(L)  •  KOGI L )  ♦  1 

IF  (KOG(L)-LEV(L ) )  6.5.3 

5  KOG(L)  »  0 

6  RETURN 
END 

SIBFT C  ANOVA  L 1ST .REP. DECK 
SUBROUTINE  ANOVA 

COMMON  P(9.9.9).INDPX(15) .KOGI 15). LEV( 15) , LUMP! 13 >. 0(21000), LIN (15 
1 ) .VECTOR (9 ) .HEAD ( 10 ) iHOLFR ( 22 ) .LEVEL  .M.N »MF 
DIMENSION  SS ( 1 5 )  .SSS ( 15) 

A1  FORMAT  (1H130X.1 0 A6 / 2 OHO SOURCE  OF  VARI AT ION6X « 18HDEGREES  OF  FREEDO 
1M7X.1AHSUM  OF  SOUARES20X.7HF  RATIO/1H  ) 

A2  FORMAT  (1H  ) 

A3  FORMAT  (8X.2A3.I22.E29.8.F26.3) 

AA  FORMAT  ( 1H06X .8HRESIDUAL 1 2 1 » E29 .8  ) 

A5  FORMAT  (7X.5HTOTALI2A.E29, 8) 

A6  FORMAT  (1H  I2.3AH  WAY  INTERACTION  SUMS  OF  SQUARE 5  -E17.8) 

A 7  FORMAT  (1M  I2.A5H  WAY  AND  HIGHER  INTERACTION  SUMS  OF  SQUARES  »E17. 
18) 

A8  FORMAT  ( 1H1 ) 

63  FORMAT  (1HO/1H06X.30HRESI DUAL  STANDARD  DEVIATION  *  E15.8) 

6A  FORMAT  (8X.2A3.A1.2A3.I13.E29.8.F26.3) 

WRITE  (6.41 ) (HEAD! J) ,J»1.10) 

LV  *  0 
DO  1  I « 1 .15 
SSS( 11*0,0 
SS( I )*0,0 

1  KOG ( I ) *0 
ERSS-0.0 
ERDF=0,0 
TSS*0»0 

DO  6  1*2  .N 
CALL  COUNT! 16) 

K*0 

DO  3  J*1,M 

IF  ( KOG ( J) )  3.3,2 

2  K*K+1 

3  CONTINUE 
SS(K)*SS(K )+D(  I) 

TSS=TSS40( I) 

IF  ( K-2 )  A, 31,5 
31  IF  (XOG(l))  A  ♦  A  .  3 


50 


NAVWEPS  REPORT  9001 


4  DO  7  J»2,M 

IP  (*OG(J)-3)  7,7,21 
21  IP  (LIN(J))  7*7*9 
7  CONTINUF 
60  TO  6 

9  £RSS-FRSS-fO(  1 1 
FROF-FRDFU.O 
6  CONUNUF 

IP  (MF-1)  27*26.27 

26  K-l 

60  TO  29 

27  K  -  2 

29  00  12  KA-r.M 
00  8  1-1,19 
8  IOC  111*0 
IN-1 

OP  ■  0.0 
ALP  -  0*0 
00  9  1-2. N 
CALL  COUNT (16) 

IP  (KOG(KA)-IN)  9,10,12 
1°  IP  (LIN(KA))  32*32*33 

32  Of  *  OF  +  1,0 
ALP  -  ALP  ♦  0(1) 

GO  TO  34 

33  OF  -  1,0 
ALP  -  0(1) 

IN  •  IN  ♦  19 

36  RATIO  a  (FRDF*ALP)/(FRSS*DF) 

J  -  OF 

^ 1 )HOLPR (K A ) * HOLER ( IN) *J, ALP, RATIO 

IP  I  IN-11 )  34,12,12 

34  IN  m  IN  +1 

IF  ILEV(ICA)-IN)  12,39*11 
39  IF  ( L IN (KA ) )  37*37,12 

37  IN  «  11 
IN  -  22 
GO  TO  36 

11  IF  aiN<KA))  9*9,38 

38  IF  (IN- 3)9, 9. 12 
9  CONTINUF 

12  CONTINUF 
K-M— 1 

00  18  rA-2,K 
KB-KA+1 
00  18  KC-KB.M 
IM-1 
99  IN-1 

DO  13  1*1,19 

13  K06  ( I )  -0 
OF  «  0,0 
ALP  -  0.0 
00  14  1-2, N 
CALL  COUNT  (16) 

IP  »*OG(KA)-!N)  14.19.18 
19  IF  (KOG(KC)-IW)  14.16,18 
16  IF  (LIN(KA))  30,30,28 

28  IF(LIN<KC))  98.98.29 

29  OF  «  1,0 
ALP  •  DU) 

INN  ■  IN  ♦  19 

INN  -  IN  ♦  19 

91  ASSIGN  96  TO  IER 
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50  RATIO  ■  (ERDF4ALP)/(ERSS*DF) 

IF  (  55-LV)  61 *61 *62 

61  WRITE  (6,41)(HEAD<J).J*1.10) 

LV  *  0 

62  J  «  DF 

WRITE  ( 6*64 )HOLER(KA) .HOLER ( INN) «HOLER( 19) .HOLER < KC) .HOLER ( I  MM) .J. 
1ALP.RATIO 
LV  *  LV+1 

60  TO  IER» I 5A»56  *60 ) 

58  OF  a  OF  4  1.0 
ALP  *  ALP  4  OH) 

IM  «  IM  ♦  1 

IF  (LEV(((C)-!M)60. 59*14 

59  IM  *  11 

I  MM  *  22 
INN  *  IN  4  15 
ASSIGN  60  TO  IER 
GO  TO  50 

60  DF  *  0.0 
ALP  *  0.0 
GO  TO  17 

50  DF  *  OF  4  1.0 
ALP  »  ALP  +  0(1) 

IF  (LIN(KO)  39.39.52 

39  IM  a  IM  4  1 

IF  (LEV f  KC )  -  IM)  AO. 40, 14 

40  IM  »  1 

IN  a  IN  4  1 

IF  (LEV(KA)  -  IN)  18.49.14 
49  IM  *  11 
IN  «  11 
INN  r  22 
TMW  *  22 
GO  TO  51 

52  IN  a  IN  4  1 

IF  (LEV(KA)-IN)  53,53.14 

53  IN  a  11 
INN  a  22 

IMM  a  IM  4  15 

ASSIGN  54  TO  IER 
GO  TO  50 

54  IM  a  IM  4  1 

IF  (LEV(KC)  -  IM)  18.18,57 
57  IF  ( I M-3 )  55.55.18 
56  IMalM+l 

IF  ( LEVI KC )  —  I M ) l 7 • 1 7 » 20 

17  I Na I N4l 
I  Ma  1 

IF(LEVOCA)-IN)  18.18.19 

19  IF  (TN-3)  14,14,18 

20  IF  (IM-3)  14.14.17 
14  CONTINUE 

18  CONTINUE 
lERaEROF 

WRITF  (6,44) IFR.FRSS 
N  a  N-l 

WRITF  ( 6 .45 ) N , TSS 
DF  a  SORT(FRSS/ERDF) 

WRITE  ( 6 ,63 )DF 
DO  22  1*1. M 
00  22  J-I.M 

22  SSS( I )=SSS( I)4SS( J) 

WRITE  (6.48) 

00  23  1*1. M 

23  WRITF  (6,A6)I»SS(I) 

WRITE  (6.42) 


52 


NAVWEPS  REPORT  9001 


00  24  I>1«M 

24  WRITE  (6*47)1 *SSS( t ) 

RETURN 

ENO 

SfRFTC  MEANS  LIST  .REF.DGCX 
SUBROUTINE  MEANS 

COMMON  P (9.9* 9) .INDEX (15) .XOG ( 1 5 ) *LEV( 15 ) .LUMP 1 1 5  I .0( 21 000 ) .LIP( 15 
1 ) .VECTOR (9) .HEAD (10) .HOLER <221  . LEVEL *M.N .ME 

11  FORMAT  (19H0  GRAND  MEAN  »  E13.6I 

12  FORMAT  (1X.A3.11X.5E20,6/15X#4E20,6> 

13  FORMAT  (1H  ) 

19  FORMAT  (1H1) 

21  FORMAT  (lX.A3.Al.tl .A2.A3.4X.5E20.6/15X.4J20.6) 

VECTOR  III  ■  0.0 
DO  8  X»1,N 

8  VECTOR(l)  ■  VECTORm+Om 
ELNUM  *  N 

VECTOR m  «  VECTOR m/FLNUM 
WRITE  (6*11)VECT0R(1I 
LEVEL  -  1 
LV  «  0 

DO  14  XF*ME.M 
WRITE  (8.13) 

LV  .  LV-fl 
LL  «  LEV(XF) 

DO  1  I-l.LL 

1  VECTOR ( 1 5  ■  0.0 
DO  2  X  =  !  *N 

J  •  XOG(XF)+l 

VECTORS)  *  VECTOR  ( J )  +0  ( X  ) 

CALL  COUNT  116) 

2  CONTINUE 
ELNUM  »  N/LL 
DO  3  J«1 »LL 

3  VECTOR! J)  -  VECTOR? JJ/FLNUM 

WRITE  ( 8. 12 IMOLEP ( XF )  * (VECTOR ( J ) * J«1 »LL ) 

LV  »  LV+l+fU/5) 

IF  ( XF-1 )  14.14.5 

5  XG  »  XF+1 

IE  (M-XG)  10.18.16 

16  DO  4  XM*XG.M 
IN  a  0 

WRITE  (6.13) 

LV  >  LV*1 
LF  .  LEV(XH) 

DO  4  XO-l.LL 
IN  ■  IN*1 
XS  *  XO-1 
DO  6  !»1»LF 

6  VECTOR!!)  «  0.0 
I  ■  0 

DO  IT  X*1.N 
IF  (XOG(XF)-XS)  7.15.7 
15  J  -  XOG ( XH)+ 1 

VECTOR  (J)  *  VECTOR!  J)-fO(X) 

I  ■  I+l 

7  CALL  COUNT (16) 

17  CONTINUE 

IE  (LV-50)  18.20.20 
20  WRITE  (6.19) 

LV  *  0 

18  ELNUM  «  I/LF 
DO  9  I-l.LF 
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*  VECTOR (1)  «  VECTOR)  I) /FLNUM 

WRITE  < 6 «?1 JHOLEROCE) .HOLER (20) .IN. HOLER (21 ) .HOLER <KH) .{VECTOR* I 
1I-1.LE) 


*  LV  ■  LV+l-MLF/5) 
14  CONTINUE 


10  RETURN 


) 


ENO 
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Appendix  D 


and  the  matrices  become 


%v hi (  h.  upon  compiitine  the  second  columns,  become 


I  he  third  column  is  formed  h\  the  vector  equation 
'1.1  \.t  * 
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where 


and 


-9.1  •  *.3 

13=  9.1  *9.1 
“9.2  *  *. 3 


C“  =  9.2  •  9.2 


thus 


9.2  *9.2 


a 

i=  1 


\2  .  2  2  2 

9„|  *(-»  *0  *1  -2 


9.1-  *.3* 


3 

2> 


il*ij=lxl  +  1  x  4  +  1 


i=l 

3 


9.2  ‘  *  3- 


-14 
C13=  3 


2912-13-  1x1  +  0x4  V 

«'=  1 


c  =~8  = -  4 
23  2 

and  the  matrices  become 


’1 

-  1 

1  - 

‘I*1 

-  4  (-1)* 

1 

0 

4  - 

¥«* 

-4x0 

.1 

1 

9  - 

¥■> 

-4x1. 

"l 

-2 

0- 

I"1 

-  4(-2f 

0 

1 

0  - 

11  x  0 
3 

-4x1 

.0 

0 

1  - 

o 

X 

Zb 

-4x0. 

upon 

computing 

the  third  columns, 

’l 

-1  i] 

l  - 

0  10' 
2  I 

1 

°  -1 

0 

1  -  4 

1 

1 

0 

0  1 

This  completes  the  Gram-Schmidt  construction. 


x  9  =  14 


1x9  =  8 
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I'hr  left  matrix  will  be  normalized  to  obtain  the  P  matrix,  and  the  right  Matrix  will  become 
the  l>  matrix.  The  first  column  is  formed  bv  the  equation 


and  the  matrices  become: 


and  the  matrices  become 
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Appendix  E 

TENSOR  PRODUCT 

The  following  equation  will  be  used  for  illustrating  the  product  of  three  tensors 

4  M  -e,, 

where 

1  <  »  <  4 
1  <  j  <  3 
1  <  a  <  4 
1  <  t  <  3 

The  elements  of  A*  may  be  arranged  as 


A  °3 


°1  °3 


L°*  “3 

The  elements  of  may  be  arranged  as 

*1 

M[=  b l  b l  bl 

bl  h\ 

m  - 

The  elements  of  £iy.  may  be  arranged  as 


11 

C12 

C1S 

21 

C22 

C23 

31 

C32 

css 

41 

c4J 

C4S 

59 
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The  tensor  Qft  may  be  formed  from  the  product^  where 


- 

a!  *n 

+ 

a? 

*12 

+ 

al 

*13 

+ 

ai  *14 

',2 

Mi 

a!  e2l 

+ 

a? 

*22 

+ 

o? 

*23 

+ 

°?  *24 

*it 

m 

a!  *si 

+ 

a? 

*32 

+ 

a? 

*33 

+ 

at  *34 

*2X 

m 

a\  *n 

+ 

°2 

*12 

+ 

o* 

*13 

+ 

«2  *14 

*22 

m 

°i  *21 

+ 

of 

*22 

+ 

a  ^ 

*23 

+ 

a2  *24 

*22 

m 

*31 

+ 

al 

*32 

+ 

al 

*33 

♦ 

a\  *34 

*3X 

m 

al  *11 

+ 

al 

*12 

♦ 

a? 

*13 

+ 

aJ  *14 

*22 

m 

*21 

+ 

al 

*22 

+ 

a? 

*23 

+ 

a3  *24 

*23 

- 

a\  *31 

+ 

al 

*32 

+ 

a? 

*33 

+ 

°*3  *34 

*41 

3* 

-i  *n 

+ 

al 

*12 

+ 

aJ 

*13 

+ 

al  *14 

*42 

ar 

ai  *21 

+ 

al 

*22 

+ 

aJ 

*23 

+ 

al  *24 

*43 

m 

ai  *31 

+ 

al 

*32 

+ 

aJ 

*33 

+ 

a\  *34 

Instead  of  having  formed  the  tensor  product  Q(j  first,  the  product  £  could  have  been 
formed  and  then  the  tensor  D}t. 

A  third  way  of  arriving  at  the  tensor  JJst  i«  to  perform  first  the  tensor  product  J,  ft'  .  This 
tensor  product  has,  for  the  example,  144  elements,  some  of  which  follow. 
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4‘fli  ^ 


•i*i 

•5*1 

•1*1 

-1 *1 

•!*? 

«1  *5 

-1  *5 

-1*1 

•!*5 

•5*5 

-1*1 

-1*1 

°\bl 

•5*5 

*1 

-1*1 

a*  62 

-1*1 

•1*5 

°!  *>l 

-1*1 

-1  *5 

•1*5 

•}  bl 

-1*1 

-1  *1 

-1  *1 

°{  ^3 

«1  *1 

-1*1 

•!  bl 

«i  *1 

*1 

•1  *5 

•5*i 

-l*i 

-1*! 

•1*! 

•5  *5 

-1*1 

-1*1 

•1*5 

•5*5 

-1  *1 

°3  *1 

•5  *1 

•5  *i 

a2  *1 

'•5  *1 

•1*! 

«5  bi 

•  •  • 

*1 

•5 *1 

°l  *5 

a*b\ 

•5*1 

•34*! 

•:*! 

"5*5 

•5*5 

*5 

“1  *1 

Multiplying  the  above  tensor  by  Q(.,  the  tensor  Qst  is  obtained,  the  elements  of  which 


‘11 


»1  M  cil  +  °*  b 2 


+  oj  6Jc12  +  a}6’c,3  + 


2  b\  en.  +  a: o; c 


°I  6i e 21 


’1  "1  t22 


2  l2  ,  „!jl3, 
ti,  C„  +  °1  61  c23  + 


,l*»e31  +  °l*l C32  +  0 


i{6{c4,  +  a 


3  i3 
1  61 


l*?c42  +  a\b\c 


C  33  + 


43 


*  -i 


are 
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di2  -  a\blcll  +  a\b\CX2  +  “1^*13  + 
aji‘c2,  +  aj6|c2J  ♦  «*6»c„  + 

°1  *2  C3l  +  *1  ^32  *  a?45C33  + 

-I**C4»  +  al*2C42  *  °J62C43 

d\i  “  a!*]cll  +  a!*3Cl2  +  a!*3-|3  + 
°1  &3  ^21  +  al*3C22  +  °1*5C23  + 
°l*3C3l  +  °163c32  +  fll  bl  c33  + 

°1*1C41  +  al63C42  *  0 1  *3  -43 

J2l  »  °2*}cn  +  a  2  '  1  -  1 2  +  a2*l-13  + 

*1  *1  C21  +  a  2  bl  -22  +  a2*?C23  + 

fl2*|-31  +  a2*lc32  +  a2*lc33  + 

a26lC4l  +  a*2bic*2  +  a$*l-43 

d22  =  a2*2cU  +  a2b2Cl2  +  "  2  *1  C 1 3  + 

°2*2C2l  *  a2*2'22  +  a2*2c23  + 

a|6jc31  +  °2  C32  +  °2*lC,t3  + 

a2*2C4l  +  a2  *2  c42  +  °262C43 


d23  "  a}*!-Jl 

+ 

°  2 

*  3  -  12 

+ 

°2  *3 -13 

a2  *  3 C21 

+ 

*> 

0  2 

b  3  c22 

+ 

°2  *3  c2 

°2  *3  C31 

+ 

a5 

*3^32 

+ 

„  3  ■  3  „ 

0  2  *  3  C3  3 

°2  *3  C41 

+ 

*3-42 

+ 

fl4i3 
°2  *  3 -43 

1 


i 


strrv*'?*  v,~ 


NAVWEPS  REPORT  9001 


64 


NAVWEPS  REPORT  9001 


ai  6lcn 
al  bl  c 2\ 


al  blcu 


a*A  bl  C41 


+  ai  *3  C12  +  °i  63  C13 


*  fl4fc3C22  +  alblC23 


+  al  blC22  +  a4  bl  C33 


*  b*  c  .  „  +  al  4?  c  . 


+  a463C42 


*4  3  -43 


+ 


+ 


+ 


A  numerical  example  is  given  to  illustrate  the  above  operations  of  tensor  product.  The 
elements  of  each  tensor  are  arranged  in  matrix  form.  The  numerical  values  of  the  elements  are 

given  for  the  Bj,  and  £ (..  tensors. 


A 


i 

s 


SI 


c.. 


~‘i 


1 

-3 

-1 

2 

aTS 

2 

TTl 

l_ 

-l 

-1 

3 

2 

2vr5 

2 

2  V~5 

1 

1  - 

-3 

2 

2  Vs 

2 

2  Vs 

l 

3 

I 

1 

2 

TH 

2 

iTs 

1 

-i 

1 

it 

V2 

V6 

1 

-2 

0 

1 

1 

1 

V3 

T\ 2 

16 

1  2  4 

3  3  4 

4  5  7 

6  7  8 


As  has  been  outlined  above,  there  are  three  different  ways  of  arriving  at  the  Q s(  tensor. 
The  first  method  will  be  given,  where  the  E ti  tensor  is  formed  first,  and  then  the  tensor  is 
computed. 
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The  Eti  tensor  is 


7  3 


Vs 

10 

1 

V3 

y/2 

16 

3 

yf 3 

21 

2 

/a 

V2 

l 

V6 

1 

sfe 

i 

V6 

0 


The  formulas  for  the  elements  of  this  tensor  were  given  previously. 
The  QJt  tensor  is 


27 

V3 

9 

2  V2 

3 

2  v^6 

24 

Vl5 

_^1 _ 

2  ’/TO 

-3 

2  VlO 

1 

V3 

1 

iSi 

-1 

2  V6 

-2 

Vis 

-7 

2  VTo 

-1 

2  -/30 

TTie  formulas  for  the  elements  of  this  tensor  were  given  previously. 
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